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ABSTRACT 


The properties of eleven types of matrix 
sequences are studied towards the development of 
fast matrix inversion algorithms. Among those 
sequences are sequences of special Vandermonde 
matrices, row-wise Kronecker matrices, composite 
Matrices, etc. 

An m-plex Kronecker system of equations is 
Ssudiedre Aprastpsolution of therm—-pilex aon turns 
OUP tosberamusualemultiplication offa rectangnlarized 
coustantevectorewith,.Kronecker®factor: matrices inver- 
ted and transposed individually. 

This fast solution of a system of equations, 
together with the fast inversion of row-wise Kronecker 
Matrices and special Vandermonde matrices, provides a 
technique of economically solving a system of thousands 
Oimeduatlonsmarisingawith the optical inverse scatter— 


ing problem. 
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(eieerPODE Colminversce woCattering Froplem 


A coherent laser beam passes through a semi- 
transparent, weakly scattering object and forms an 
interference pattern under far-field conditions. 
ismiieepossible te determine the optical structure of 
SUCH waNneObjeCcuemMaking Use) of the intormation recorded? 
This is the optical inverse scattering problem. 

An interference pattern records only the magni- 
tude of an optical wave. The phase of such a wave 
cannot be recorded directly by any instruments of today's 
technology.) Without thas information, the spatial dis— 
tribution of the phase of a propagated wavefront, and 
Nencemtnewoparcalestructure Of (a Scatterer, Cannot be 
LeconseLuctcd mr Orethisereason,s Lie opeical anverse 
scattering problem has been considered a difficult if 
(1OtronmUnscoLVablesone:. 

A breakthrough was made by Gabor (1951) who re- 
constructed the image of a three-dimensional reflecting 
object from a two-dimensional photographic recording. 
His technique was later developed to become holography, 
and convinced people that some phase information was 
indeed recorded in an interference pattern. Contribu- 
tions Gue to Goodman (1967), Lesem, Hirsch, Jordan 


(1967 781968), Hickling (1968) proved that a three- 
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dimensional image can be reconstructed not only by 
optical methods but also by computers from holographic 
data. WoLe (197 0l—=1 97 )) anc Carter (1970) studied 
three-dimensional transmitting objects by off-axis 
holography. Wolf concluded that some of the spatial 
frequencies of such an object can be detected from 
nologvaphie@eatasbutwthat an=order to achieve a com- 
plete three-dimenstonal reconstruction, thevobjecteunder 
study has to berrotated. ~He pointed out that ‘hologra- 
phic reconstruction may be good to a resolution of about 
9 times the vacuum wavelength of light. 

Another eae ch towards the inverse scattering 
problem parallel to holography went on in x-ray crystal- 
lography. X-ray radiation has a much shorter wavelength 
than an* optical laser. ~Kendrew (1962) “and Perutz (1964) 
attached a known molecule to the molecule of unknown 
structure, to measure both magnitude and phase of the 
diffracted x-ray radiation, and were able to sketch the 
profile of the object molecular structure, by comparing 
the interference patterns obtained with and without the 
known molecule attached to the object molecule. Sayre 
(1974) claimed that computer reconstruction from x~ray 
patterns carried out by IBM researchers achieved a 
spatiale resolution or 1.15 APemBoEn holography and x-ray 


crystallography represent a diffracted field by uniform 
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plane waves and process both the magnitude and the phase 
transmitted by these waves. 

There is a recent attempt to measure and to pro- 
cess the magnitude and the phase of plane-wave compon- 
ents of the nonuniform type of the field arising from 
a diffracting or refracting object. This new approach 
May be divided into three steps: 

A. To recover from interference fringe patterns 
the® spatealy d2stribution or magnitudes and 
phase of an optical-wave field in the far- 
field region. 

Boe lo compute Enos spatial@disteributcion of the 
object near field from the complex-amplitude 
field data in the far-field region. 

Gye Tomdetermine. the optical Structure of a 
Scotter ngeobjvecterrome ete. near-rtieldrdata. 

~tep A was solved’ by Schmidt-Weinmar (1973(1), 
1975(1)), who proposed a three-reference-beam method. 

By this method, the magnitude and the phase at any point 
in the far-field region can be computed from three 
interference fringe patterns, obtained with one and the 
same object field superimposed upon three versions of 

a coherent reference field. Schmidt-Weinmar (1975(2)) 
abso published agtheory concerning step B. Step © has 
noc yet been completed, but preliminary results concern- 


ing step C were presented by Schmidt-Weinmar (1971, 
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2273 (eto nesOprca laSocitety, of America, 

With the new approach, the scattered field at any 
point in the far-field region is considered to be the 
Glove leet toctsor manyehuyden's point source fields, 
that emerge from inside the scattering object. These 
SOuULrcemiicldc. beingractivated by the ancident laser 
beam, may number in the thousands. In mathematical terms, 
Phewspetialwecistriburrzonvworethe, scattered far -fieldis 
petabcds@ a .tlescpattalidicstributwonsok these point 
sources by a linear mapping. To determine the spatial 
distribution of these point sources requires the solu- 
tion Of a linear system of possibly thousands of 
equations. This system of equations is neither sparse 
Nor diagonally dominant. Solving 21t by Gaussian 
Elimination is impracticable owing to the hundreds of 
pages of storage space and the amount of CPU time 
required. Hence, a special computing technique 
is sought that requires less storage and computing time 
than any of the conventional methods for the solution 


of such a system of equations. 
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1.2 Conventional. Methods of Solving a Large System of 


Equations 
1.2a Solving a System by Matrix Inversion 
Se ee ee ee VS Be UOT 


Consider a linear algebraic system over the com- 


plex field 
Mg SY (Cise2 aan 


where M is a non-singular square matrix of size (oo) 
eC ei el Secateurs eo teste s(n, mM) COlizat na nO am eSe smo 


data. The system will yield 


xX = M xo Gey) 


where X has a dimension (n,m) giving the m sets of unique 
So Ucw iomcom Gi. ll) nect Moissunitary, «1.202) wild 
simply become 


= WE x Gh 


Me being the Hermitian transpose of matrix M. Otherwise, 
mt needs to be evaluated and the lengthy calculations 
normally require ‘the use of én Bieerr anc Computer: 
Apart from obtaining the solution X by inverting 
PUemCOer Er ClentLamatrix: Mya (ls221)) can bewsolved by other 
methods as well. These methods usually fall into two 
categories: direct methods and iterative methods. Each 
method has its advantages and disadvantages and selec- 
tion of these methods depends largely upon the structure 


of the coefficient matrix M. 
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ecb Directs Methods 

The major difficulty in the use of direct methods 
lies in the number of arithmetic Operations required. 
Table 1 shows the number of operations needed by 
various methods, Cor Elem Ceasccsinewit cis 1sean (nn, 1) 
vector, as summarized by Westlake (1968). With the 
exception of triangular and tridiagonal systems, all 
direct methods need a number of multiplications propor- 
ELOnAal “to i WHEN Ves large, say, mis darger tnan 20, 
direct methods will be costly. The problem will become 
particularly serious when M contains many elements whose 
absolute values are less than one. Repeated multiplica- 
tions of these small elements lay sresultein underflow. 

Another weak point of the direct methods is that 
they usually involve too many SUDEGaACcoLolG,  wittche willl 
easily magnify percentage errors. Nowmaiad ya darect 
methods are considered applicable to small systems of 
equations. For large systems, direct methods are used 


only when no other better methods can be found. 


1.2c Iterative Methods 
OnewoL Une intrinsic advantages of these methods 
is that errors due to computer round-off etc, may be 
damped out as the iteration continues. The major dis- 
advantage of all iterative methods is that convergence 


usually requires a diagonally dominant coefficient 
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TABLE ~-1 
Number of arithmetic operations 
(A is an mn Xn matrix) 
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Method aD 
Direct = 52 + & = 
Gaussian 
Elimination n tn? +n? — In dn? + nt —2n 
Jordan 1 tr + n? — bn tn’ — in 
Doolittle n qn’ + n? —4n $n3 +- ni? — $n 
Cholesky n an? + 3n? + dn an +n? — in 
(symmetric) root 
. recip 
Crout yi? + 32° + gn 
Tnangular n in? — tn $n? — in 
system 
Product form tn’ + $7? 
of the inverse 
Tridiagonal 2n 4(n — 1) 4(n — 1) 


— rr ne ae ee se 

















Jterative Per Iteration 
eee eer eer eee Sree eee 
+ x + & = 

Jacobi and - 

Seidel n n? me? — 7 
Steepest 

descent 

(gradient) 1 2n? + 3n 2n? 4+ 2n —-2 
Conjugate 

gradient 

(not symm.) 1 3n? + 6n +2 377 + 32 — 3 
Successive : 

overrelaxation n n(n + 1) ni? 
Peaceman— Sn 

Rachford n 2n*? + 8n — 8 2n? + 82 — 8 
Newton- ; 

Raphson 

Analogue 0 2n° 2 2h art 





Matrisc. wer Ornra asuiitab le vclass of Systems of equations, 
iterative methods are very efficient, while for other 
classes, iterative methods need not even converge. 
When Y is an (n,m) matrix, both direct methods 
and iterative methods have to be applied m times in 
ordeveto Getmihem sets tof answers forsx. Therefore, 
the procedure of inverting the coefficient matrix 
preliminary to solving the system of equations is 
often favoured, since we only need to invert M once 


fore al in seroeoLedata. 


1.3 Redundancy Techniques in Matrix Decomposition 
Se ee ee i ee Se SPOS ye On. 


Good (1958) showed that we may encounter matrices 
that can be decomposed into a product of several sparse 
matrices, and that in this case manipulating the sparse 
matrices can save a number of mathematical operations 
and storage space. He suggested the technique can be 
applied to speed up matrix-vector multiplication. 

Cooley and Tukey (1965) aetiner ster created 
redundant factors in a Fourier transform matrix so that 
the matrix can be decomposed into two or more sparse 
matrices for easy manipulation. This is the well-known 
Fast Fourier Transform algorithm, which paved the way for 
déVeElLopmentnot VeteqgvastsHadamandeirancsctorn pnetonpedGea 


later date. Theilheimer (1969) formulates the actual 
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procedure for explicitly carrying out such a matrix 
decomposition. 

Andrews and Kane (1970) applied the techniques 
to build a Kronecker matrix model, and they proved that 
many fast transform matrices are in fact special cases 
of their Kronecker model. Their investigation has led 
to the discovery of new orthogonal matrices such as 
generalized Hadamard, generalized Walsh, and generalized 
Haar transforms. 

The advantage of employing redundant factors was 
summarized by Andrews (1970): suppose that a matrix 
of order p" by ane can be expressed as a product of n 
sparse matrices each of p- non-redundant entries. When 
Bne*Orloinai=-mavrix is mureienied by -a VECtOT "on size 
po only np tt multiplications and the same number of 
additions are necessary. This compares favourably with 
the pee multiplications and additions normally required 


if the matrix is not decomposed. 


1.4 Solution of Inverse Scattering Problem by Matrix 


Decomposition 


tHessystem of equations that arises from Scatter-— 
ingethicory (see Eq. (7.3./)) may be put) in the form 
P!'f£= F', The vector f is the unknown, which contains 


as its elements the Scattered-field source densities at 
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points (XYZ) Vere, 27% «0, Pwithini tthe scattering 
object. The vector F' contains the scattered-field comple x 


amplitudes at some points (K_, 


a eras Ie eee. ae 


Ryu Keu 
inethe far-field region.s! The Prepagalor, MatrixeP ‘Macts 
as a linear operator mapping the vector f into the vec- 
comer ee Withiia givens sety of data F" and sas known propa- 
gator P', it is possible to determine f by solving the 
system of jequations P')f = F'. The present problem is 
how to obtain f quickly and accurately. 

Both P' and F' depend on the two independent 
variables Ky and ye whereas £ does not. We are free 
to choose in the (KK) -plane SOmMespanreVen bai pointes 


(K Pawar ciawi ie giverPiea Simplatied form’ P mand 


rk 
xu Vu 
modify) F Waccerdinglyato F) to-aldowla fast solutionsof 
Pa eee nem tara Sldedubys Patu—al iwi Il be the same as 
ise) Ge Nakeaveia. tops: 2) oe egg 


P' has all its elements expressible in terms of 


the products of three periodic functions: PS = 


exp (-ix. K Nee ya ry Hb ene eae If we make Se 


View 
take the same values modulo 27 for all elements in any 
column of P, then we have created a redundant factor 
exp (-ix K,). Lita Slider smannie nye Lt YK yu has the 
Same value modulo 27 in any column of P, we have a 


redundant factor CoD eae et) These facts enable us 


to apply decomposition techniques to handle P. 
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Pogeactuaie practice, let us assume that our 
optical problem requires an equi-spaced grid of source 


points. We can write ce) ere nee ings Sal Boo) 33,42, 


Wietvesx hy. 42 are constants and Iye Jor 33 are three 


Niel 


SOs Ole positivesintegers., | In order to sample F regularly 


ate points ion OS plane, we shall assign aes = 


KAR, k,, (2p7/Ax) and Sih = a + ky, (4am/Ay), 


where k k k Kar Dr mo ewe) aeintegers and AK, and 


Tey OY Cia 
SO aCe COns tates. By letting ae Jor 33 be increased 


ROWoWESom lox LCograpnicatlyein steps of )1), and k k k 


2 aM, 


k be increased column-wise lexicographically in steps 


“al Uf 


ofel;, P will become a-matrix of row-wise Kronecker pro- 


ducts of three smaller matrices R, S and a with m= 


EZ feck NOL e= kh, So, eand a correspond to the three space-~ 


sented functions in the elements of P. 

cl ome UMC Cuca WOU DoInted sour that Rk, Sand 
Tn all are of the Vandermonde type and we already have 
quick methods for inverting them. Hence we can expect 
a solution of the system of equations P f = F results 


A oil 


ciethe seOcimece—— (oaGuR \ um @ ar ype. where F" is 


aeVECLOG LunctLion wor Nae 
ingtnhescourse of developing oursfast ysolution for 
Plt) = lj.a Nunber of fast algorithms weresdiscovered 


for the inversion of special types of large matrices. 


We have found that among the Vandermonde matrices, there 
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are two special types whose inversion requires only 
conjugation or physical rotation of the matrix. It 
is convenient to make R and S of these types to save 
arithmetic operations. 

We have further noted that the Kronecker CDECOGUCE 


vi 8 RS) is still a matrix of considerable size so 


that matrix-vector multiplication (ca ® Ro+)pn becomes 


(S 


lengthy. Improvement is made by displacing side-wise 
the elements of vectors f and F" to form rectangularized 
arrays (Rect f) and (Rect F") respectively. Then system 
Evre- lis etinally transtormed into Rect (Rect £) = 

Re. (Rect EY) (ees \o The unknown complex scattered~ 
field source esac t are obtained in the elements of 
the array Rect tRect £). 

This thesis will be presented in three parts. 

Part one describes the properties of eleven types 
of matrix sequences which appear in our main algorithm, 
and develops some fast algorithms for computer inversion. 

Part two deals with an application of the fast 
algorithms to the optical inverse scattering problem. 

We develop a special system of equations for the inverse 
scattering theory, and a complete method of solution. 

Part three consistsrof an exrom analysis of the 
inversion of those sequences introduced in the thesis. 


ANSerror bound us found for the methou of solution 


developed in Part two. 
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PART ONE 


FAST INVERSION OF MATRICES 
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If. SELF-INVERTING SEQUENCES OF MATRICES 


2.1 Matrix Sequences 


Definition The matrices eg an cts MM ay see 7 
mM, are called a sequence if they are related by a re- 
CuLrsion formula 


ee ee Meee ee Mo | Man) be (rsdlyal) 


O@SSpecialginterest to,us. iS) the) simplesti case of 
the above, in which a matrix in the sequence depends 


oniyeonieheematraxt just. before it,snamely, 


rg ISOS (2a 2) 


This type of sequence will be the main Subject .ofedis— 
CUSS Hone neLhisethesis. 

Definition If all the member matrices My iMoreeeee 
M. are orthogonal matrices, the sequence of matrices 
will be called an orthogonal sequence of matrices. 

Definition The inverse of a sequence is defined 
as the sequence of inverse matrices. 

Theorem 2.1 The inverse of an orthogonal se- 
quence of symmetrical matrices is just the sequence 
Ltselty 

Proof: .Since the,sequence consists only of ortho- 


gonal symmetrical matrices, which are self-inverting, 
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inverting the sequence will not change anyeoteite 
member matrices. 

me GeneralazecaDerination Sinathasethesis, a 
property modifier will be attached to aysequence if 
the modifier describes the common properties of its 
Member matrices. 

An example is the term "orthogonal sequence" 
defined above. This term means "a sequence of ortho- 
gonal matrices", Other terms to appear later in this 
thesis willbe tinvescedl, “symmetrical, "Quasi- 
orthogonal", "product of" sequences, etc., which will 
respectively contain the meaning of sequences of 
pitvetlcd ae cymmetrical™, “quasi-orthogonal". PDEOUuCE 
of" member matrices, etc. 

Melero Oe cetonceme (le Lhd suthesic ab eiss often 
necessary to break up a matrix into its components, 
OretOerotatesi tt by 20 Udegrees, or to turn) a matrix up 
side down. For the sake of convenience, the following 
notations will be adopted: 

M(n,m) represents a matrix M of dimension nxm. 

M(n) will be used in the place of M(n,n) for short. 

M(n) will indicate a row vector of n elements. 

M(n') will represent a column vector of n elements. 

Mt means the transpose of M. 
mu means the transpose of the inverse of M, i.e., 


(M71) *, 


is 
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Se 


indicates the matrix M being turned upside 


down. 


=> 


indicatecethe Matrix Mi being physically rotated 
clockwise through an angle of 90 degrees. 


For example, 


Teak) eRe alee 
M(2,3) = 
Olean ee 
eye os 
WO = 
aie, ea 
Soy Vie 
MCS aA) eS Mo 5 My 5 
Ase ake 


Whenever possible, we shall use lower-case letters 
apo; GC, eeeetO denote Glements of a matrix. These small 
Petters wai be bessubscripted or soninie aaa SNe to 
indicate the row and column where an element stands in 
SRVCCUOMIOMeIn aumatrix. lhe: Copltal’ lebbers spb hep ee. 
will be used to denote a matrix or vector, and will be 
SubScripted if they. are rows or Columns Of @ matrix. 


Thus, following the above examples, we have 


M, (3) = [m,, ™o2 ™ 31, Moa") = : 





| Loe 

tied Taal Sas eg teataMll,t = 
* 7 i 

| = + spe aed $6.0 sai oe _ 











»£4- . i f a 
, eon « \ : mf we , ae 
. ; ; 
ii 
ij 
| 
i 
sf < - Uo we Log ‘) 
EC Lena j ve 
4 : Gis 
7 py’ SD ayl es [ hay : 
T..7 


at Saban. 5 wows: iti Li oOy 
tk. (2ykae Beg tek.” inc? fag? 3S Bhos 


Aw" Brad cea ee! i 





_ ; 


a>. s 


est oF 


The double subscripts will be written simply 


Side by side, as in m to save typing Space. When- 


Zoe 
ever confusion may arise, a comma will be used as a 


Separator y,easiin m When a Variable, ¢€.9g. r, is 


A eed 
subscripted with another subscripted variable, ere) Kos 
TOE Jor Bie mSUDS crip of thersecond variable or 
variables will be placed along side for convenience 


insty pring. Thurs ry will be typed as r as 
2 


k2 35 
Y Dil, ae Sie es 

BULENemmnore, sagDresi xX subscripe attached toa 
Matrix or vector notation would indicate the label of 
aeuUunberede mati smo Vector... Hence 

yM (2,3) MeaiSmilaurl eNO. slEOresize (2,3), 

4M, (3) means the second row of matrix No. 4, 

with 3 elements in this row matrix. 
2M, (2') Mees ene thitdscolumn of matrix Nos2, 


With a size of 2) elements. 


2.2 The I' Sequence 
Definition Sequence I' is the sequence gene- 
rated by the recursion formula 
O(n!) TEN A) 
NEC a ea alley = J Pe ey Ole eee ee Lt) 
1 O(n) 
with sae (Ls areca [tlt ee 


Sequence I' is also known as the "Exchange Matrix" Sequence. 
Properties 
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member matrices I'(n) equal to the mirror image of the 


Pdentity matrices .—(n)). Thus, 
1 
O 
ih 
TE Cas ail 
le O 
(22) 
(y= te (ny) 
WEL wEeindi caving 4 transpose. 
2 oeebeeece an gor thogonal. sequence 
1 nt — f 
16 Y a Gey yak ce. = DRG Se (OE 2 )) 


See lie aDovesiidicatessthat the Square of any 


member Matrax—of I’ is an identity matrix. 
1 2 — 
rere (2) a) ieee (1) (200) 


ee cCOMel icoremu sly tats Cleat iia Sequence 


I' is self-inverting. 
ea = 1 
it Ray) Ss aly (29275) 


5. Multiplication .-by Sequence h' trom lett 
will turn a sequence up side down. Let the sequence 


be M. Then 


I'(n) M(n) = M(n) . (Seore) 
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For example: 


Opel tA ie) a) 
: OP PAD) 


iF Ke) ZS) 


6. Multiplication from right by Sequence I' will 


turn any sequence member left to right. For example, 
2a Ome Sy 
= . (rece 


From property 4 above we have obtained our 
Pas ce voor. enim: 
inewinverse of (cequence 1" iswqust the Sequence 


Pee Sele. 


2.3 The A Sequence 


Definition An A sequence is one which satisfies 
CPieereCULS on cokmiula 


(ies cele) 


A (k) O 
A(k+m) = 


O I' (m) 
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Theorem 2.3.1 Sequence A is a symmetrical sequence. 
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will also be symmetrical. Since the sequence, starts 
woth AQl)e= (1) which is symmetrical, by mathematical 
induction, all member matrices in the sequence will be 


symmetrical. 


theoreng2.o.2) A(n) is an orthogonal sequence. 

i 
Alp lAr (nye= 7 (ni) %. (Oe 30) 
PLOot ee Tome(2 2522) and) (2232 1)%) with Teme i 


A“ (k) O 
ie 2 
A (ny Agim) =A (n jie (eso) 
O r'? (m) J 


BLOM 2m 2e4)eet'o(m) = lm) which ts an identity 
matrix. Provided A(k) is orthogonal, from aS Sees 
jam tS = A(k)A™(k) = 1(k) Wiichetce also an @denerty, matrix. 


(293% 4) becomes 


tek) O 
Rcae ene = =I(kim)=I(n). (2.3.5) 
O I (m) 
Therefore A(n) isPorthogonal 2) Now A(1)) = [Ph] as ortho- 


gonal; by mathematical induction it is established that 
all matrices forming Sequence A will be orthogonal. 
Property of Sequence A Any member matrix A(n) of 
the Sequence A is a permutation matrix. Postmultiplying 
any matrix by A(n) means to permute the respective matrix 


columns. The following example illustrates the manner 
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of permutation: 
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means to permute s 
AS a special 


A sequence exists 
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case, Then the 
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for 
n= 27 Yr = a positive integer . 
The sequence is 
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A(8) = |----=--- ital etc. 
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iiesrecursion formula for this special case is 


A (n) O 
A 2 es 
O Tete (ni) 


Theorem 225 gives a 
Fast Algorithm The inverse of sequence A is the 


sequence itself. 
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IIL.  QUASI-UNITARY SEQUENCES 


3.1 Quasi-unitary Sequences 
Sd ee ne ice 


Definitions Chapter II introduced: two sequences 


I' and A which satisfy Rip aes TeandsAn = lie Seve | Wee 


and ne = I. Fast algorithms were derived for these two 
sequences. 

To generalize the algorithms, the identity matrix 
I above is now replaced by a hon-sSingular diagonal 
Matrix D. We shall call any Sequence? OF traving the 
characteristic elohe = D a quasi-orthogonal sequence, hav- 
ing the characteristic gq/= D,O = 0% a quasi-unitary 
sequence, and having the characteristic ar = Da quasi- 
self-inverting Sequence. A quasi-self-inverting sequence 
iSeaespecial™casewof A quasi-orthogonal sequence, and a 
quasi-orthogonal sequence is a Special case of a quasi- 
unitary sequence. 


Their fast inversion formulas can easily be derived 


ELom, thesderintrerons 





Type of sequence |/Characteristic]Fast inversion algorithm 





g-t=p lo, alternatively 





quasi-self- 








inverting on-1= op71 

quasi- g-l=gtp-1l 
orthogonal g7 l= p-lgt 
quasi-unitary QQH = p g-1l =QHp-1 
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Pict teem vont a diagonal matrix D, it is 
only necessary to reciprocate its diagonal elements. 
theverore, the algorithms in the above table indicate 
that to invert a sequence Q, we need only to transpose 
Q (or conjugate and-transpose Q), and modify its rows 
or columns as appropriate with the reciprocals of the 


respective elements of D. 


3.2 The Diagonal Sequence 


Delinmerony Foequence D 1s any sequence satisfying 


the recursion 


Acai) = (Gaze) 
O D (m) 

with k and m being arbitrary positive integers and D(k) 
and D(m) being diagonal matrices. 

Properties The following properties are obvious. 

1. Sequence D is a diagonal sequence. 

2. Sequence D is quasi-orthogonal. 

3. The inverse of D is sequence mes another 
diagonal sequence with 


Ages only dae ot We (see 7) 
uu uu 


Fast Algoriuthm To 1nvert D(n), 12t) is only neces- 


sary to reciprocate the diagonal elements. 
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Jeet hie lB Sequence 


Definition Sequence B is defined by the recursion 


formula 
esp) Ae) I (n) L{n) B(n) B(n) 

B(2n) = = Corr ede) 
O TCT eee cy) == (1) fn) 1) 

Cees UCN ap 


It follows from the definition that the B 


Sequence consists only of those B(n) where 


i pee ¥ = a positive integer. (oe) 
Therefore 
BCL eee) 
LOS 431 meee, ie oR 
BZ) = — 
OP se — 1 -l 
eoete Le | 
[B(2) O 1 (2) WT 25 1-1 1-21 
BC aye = (37353) 
[ O C2) Wel (2p ert ( 2) oO. — 1250 
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Pees Os Oe Ol OF 10 
Omer een Olen a Ol 15 Oo 
Op Mey ak de) 6 tek = xe) 
OR Ome Oe alee Cant: pO aad 
Sean tes 


pequence Bis clearly not orthogonal; 
MHCOLeEMes.o poeguence B is’ quasi—orthogonal. 
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B(2n) = = (S334) 
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(S43 oy meelC1 Catesr that ant B(k) B® (k) is diagonal then 


B(2k)B* (2k) will also be diagonal. Since sequence B 





Cub 


starts with a diagonal matrix B(1) = [1], it can be 
Seen from=(3.4.0) “that Bie Bs (2)ea(ajBo (4) ies. are 
all diagonal. Therefore sequence B is quasi-orthogonal. 


Now let 
Bin be (eee in) (2B ee) 


pe LUingeki—el), 2,4, oye. 6iny (3.3.5), the actual value of 


sequence D¥"1s calculated as follows: 


DCL) ieee | ds) 
Z 
Da 2) ee 
2 
4 
4 
D' (4) = 2 
2 
(Sensei. 
8 
8 
4 
4 
D'(8) = 2 
2 
2 
2 
The recursion, formula for sequence Db!) as 
DEFT G a) 8) 
DIS C21) eee ° (3,2.0) 
O ata) 


PUCSUCSRINVErSe. 1 Seround gin e(3 2.2). 







“ 5 : 
: ' 
s 


dims bY oft? 2-00 ee Dae oe 
sen tu, et ved Cg aRS Rats ee 


. A 6k Oe 6 nm ie le e as | eh nd SA oaty 4 ces 





Py » (Be 
q 


i 

| a 
| 5 L 
' ° 7 

\ } 





e? '\ sumhgpss, 20%) a lupeos ahes wos 
7a 7 A 





28 


Fast Algorithm The inverse of D' can be obtained 
by reciprocating the diagonal elements of its member 
Matrix. Or alternatively, reciprocating the first half 
of its member matrix (diagonal slementesonly) +) and 
replacing the diagonal elements in the second half with 


aeConstants) 5% 


Dimon) = (83595 
O 461 (n) 


Then the inverse of Sequence B is obtained from 


(Snes. 0) 


as 


Bee Eee EL (nD Lal) Ch aalis 


In other words, the inverse of Sequence B is the trans- 
posed Sequence B with each row of the member matrices 
modified by the respective reciprocal of the diagonal 
elements of the member matrices in sequence D'. 
Sequence B can also be called a non-normalized 


"Haar matrix" sequence, with proper permutation. 
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IV. VANDERMONDE SEQUENCES 


4.1 Vandermonde Sequence 


Definition Given a sequence of distanct numbers, 
real or complex, denoted by XyrXorXgreee se Xye a sequence 
of matrices, called the V sequence, can be generated by 


the formula 


Way) X, (n') 
V(ntl) = ol AIEEE CL aa, pero (Are) 
- n 
AGN esa ace) 
with 
VG) eee | 
-, sat key Bal n 
X, (n) = [xy Ky Xq eeeee a 
ts 
2 3 n-l 
q aut 4 
SAG VES Galt ee eer cer! 


The member matrices in the V sequence are called 


Vandermonde matrices of non-confluent type. They take 


the form 
ws i al a ar te te eel. 
xy Xo X3 oe 2 X 
a ce x4 Ha ye Fs cca 
Viti) == (ae) 
xncl ynnl ynol xnol 
1 2 3 a 9 
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Macon and Spitzbard (1958) first derived a set 
of formulas relating a Vandermonde matrix with its 
inverse. These formulas are complicated. Parker 
(1964) worked out another method SubeaDleslor computer 
implementation. On the same basis, Traub (1966) designed 
twosatcorithms.. He, recommended his Algorithm I. for 
De- ee ew ch requires n(/n—-9)/2 multiplications and 
brie l)Y ZvadcrlLions to anvert a Vandermonde matrix. 
Ballesterm and Pereyra (1967) found that to solve a_non- 
confluent Vandermonde system of equations, only n(3n-1)/2 
multiplications and the same number of additions are 
required. However, using their algorithm to invert a 
Vandermonde matrix requires more arithmetic operations 
than using, Txraub's algorithm I. More contributions on 
the inversion of generalized and confluent Vandermonde 
Metin seo bemOUem LO han LO.) ebgorcke and Bliving 41972), 
SECU D pe wow oy 2 andeGoknare (1973) . 

In this thesis we shall follow Parker's formula 
and Traub's algorithm I as we believe that these are the 
Simplest for the inversion of non-confluent Vandermonde 
matrices. We also develop algorithms for special Vander- 
monde matrices. 

Inversion Method A simple method to invert a 
non-confluent Vandermonde matrix was due to Parker 
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Let a polynomial G(x) be formed with roots x, 
(Vee yt) edual sto the elements in) the second 


Boweo fave) inet 21. 2h): 
n 
CCS) eee (a ee) Arley) 


Then let g,, (*) Demanotniewepolynomddaleol x such that 
Q,,(%) = Glx)/(x,- x) . (4.1.4) 


Compute the coefficients of ie a IO ACN C I17, Gi(x) being 
the derivative of G(x). These coefficients are the 
Clements siieLowsll Obea Matrix (1). | Proceeding with 
WS Meee Sa odepityy Gl SE ikhenenboue by Gap ls seul be formed. 
Then M(n) is the required inverse of V(n). 

PLOGt eo cece eMethiod: 

VaerOoMmene  COMpitarlOn OL Min) and 421752), it 
PseecleCaretnacerne product M(n)V(n) will be a matrix with 
its element in row u and column v equal to “J, (8) /G' (x). 

2. From the definition of a derivative, we have 


Gi) im G(x) / (x-x,,) ceed gas 52 


XOX 
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Mate GeWelete lnehets, toWeleuvepe (Wiel Sy) | om eas ) according to 
(4.1.4). Therefore G'(x) approaches -9,, (x) when xX 


approaches Xr Or 


“G(X )/G' (x) = 1 smoke Suh C2 1 (AiG ) 
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9, (x,) = 0 or Sh RAS Aes = OM OG eV eu (a) 87) 
4, even 4. 6) ands(42,7)) at follows 
that 
MA yeeViAt)) =). 1e(n )P 2, ae) 


Mat .Seto say 
Mi (ces Vile) (4.1.9) 


This gives a general formula for the inversion of 


(ATE? Zanes 


vt(n) = Din) vin) (4.1.10) 


where D(n) is a diagonal matrix equal to 
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ahs hey BAI, ) 


V(n) is a square matrix whose elements are elementary 


symmetric functions 
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Fast Algorithm Traub (1966) denoted the elements 


of an inverted Vandermonde matrix by 
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He then suggested that any jth elementary symmetric 


function can be calculated by recursion formulas 
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plications and the same number of additions. 
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All the ne may bevobtained with n(n-2) multiplications 
and n(n-1) additions. Then q.; may be computed as below 
for n(n-1) multiplications and same number of additions: 


= q Us ae, 

wie = cee T ap bh teak ives, Olle «ite 

Vel 1812 jetty (ar eieec) 
ue ' aL 

St oe Mea ey 


Let us assume that division time is comparable to 


aemuibLociplication tine for the above, calculations. Then 









i an e 7 
ou ee _ : i / 
‘ 74 vhs oo ae 7 _ 
oye -F eat hta =U jae , neL Hite 34 ss 
. 4 we +. an : ae 
an (peta? a lcAgne he vd Pee hAL its a br a ae 





7 = “ Fi Pe ad 
(Ups ; ot ae fe ati Ve er 
Fae 


' 

, 
a 
a 
e 

a 


9 r: ; ae 7. te fT ) Fes 7% 5 sagt i aC eto) is f Ny ah ae \ oe ‘ce 
we | syci Be fy es + cit 12 eh oa | " +, ea! ‘ PRG te ( Matte ae La 5 
pte Sethe Yo. 1odinen ao uo eee oe i +c urn (ks alee 
. aie. 2 h ’ 





Ciewed> ee rc pei Traub s algorithm is to perform the 


ive otelOoncm ine (ebb )., and sthise needs m(n—L) 


multiplications (the first column of V(n) is always 


eu oavectom)r, 


plications and) 5n(n—-1)/2 additions. 


Altogether, we need n({7n-9)/2 multi- 


Weighs, Shek UMEe yen eS 


first algorithm for inverting a non-confluent Vander- 


mide matrix. 


4.2 The C Sequence 
Definition The C Sequence is generated 
formula 
E(n) X, (n') 
C(nt+l1) = 
Nal 
a 
waclet eds) es [x,] 
eee Zee n 
X, (n) = [x soe eer ee ] 
es (rts) t= xo ee hie 
Z a i BA ce ee ee 


ABMengereMacl Vad ne ene 


Awe) = 


Sqo0eG5cnHO0D &X 


n n 


byiechie 


Car 2 a) 


(4.52.2) 


30 


Ror ten 
:° 
he ws 





Properties 

1. Every row in C(n) is a power sequence. 

Zee (Th) Cans letpacClOorl Zede into two mattrces, 
the first being a diagonal, one, the second being a 


transposed Vandermonde. 


2 n-l 

xy g) xy x Roe x4 
2 igre 

X5 ip Xo x5 ee Xo 

Cena) = 8 eesceceeeeee = Dinky (n) 
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xc ily oe x Aare e 
n n n n 


(Anis) 


SO wee emu VveESemOLRC (Ni etsmobtainedetrom (44253), 


namely 


= = lp Ih 


Gee aD ban) jae ly an) | MeD aM Ev ein a De (a). 


(4.2.4) 


Fast Algorithm The inverse of Sequence C is 
simply a transposed inverted Vandermonde sequence, 
Mmodiaied byVasreciuprocal Yof "thetiicst column of Zits 
Member matrix, the Vandermonde being formed by dividing 
each row of C(n) with the first element in each row. 

An example is given below. 


When n = 4, we have 
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~1/1x, (x5-%,) (x3-x,) (x,-x,) } 
“1/1x,(%4-Xy) (%37Xy) (X4-x,) I 
~1/(x, (x) -%3) (X5-x3) (x,-x,) } 
= —_ i — ‘ -_ 
1/1{x4 (%)~-X4) (X%5-Xy) (X4-x,) I 
Ces) 

Pereyra (1967) warned that attention should be paid 
to the possible instability of Vandermonde matrices as, 
generally speaking, these matrices can easily become ill- 
conditioned. In the next two sections, two special types 
of matrices will be introduced which are not ill- 
conditioned and for which there is a faster inversion 
method, These are the types of Vandermonde matrices that 


are to be used in our optical equations. 





4.3 The E Sequence 


Definition Sequence E 


Sequence C with 
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For n = 2” and r = Do aR 


the following form: 


is a special case of 
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Since Sequence E is a special case 


= exp (-in) 
= exp (-it7/2) 
xX = exp (-iT/4) 
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of Sequence C, the 


fast algorithm.described in (4.274) will be applicable. 


However, a still faster method can be developed with the 


help of Sequence I!: 
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theorem 4.3. The product of Sequence E and its 


Eranspose iseequal to the Sequence I’, multiplied by n. 
i: — q 
ALGO) Th Gay) Se aay ae tGae) (Caneee 3) 


Proor-meoccordingeto (473. |) and. (422.2), a member 


matrix of the Sequence E can be written as: 
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TiesiateixepLoduert scan be reconstructed from oF 


iy ae n = ie ED (ee (28207) 


Therefore the theorem is established. 
pecclone2.2 explained) that) £' sequence is self— 


Tivereund ss | Eromat(4.o.0)) "it 2s clear that 
-l age oa 
Hs (a) = Wee ea Aa eB (4333, 8) 


Bo (n) IsAtnj means to xnrotate the matrix E(n) clockwise 
Lorougnedserghitweangle, and no multiplication is involved. 
This leads to our 

Fast Algorithm The inverted sequence E is the 
original sequence rotated clockwise through a right 


angle@andmthien  dividedeby. nn. 


if e) & ialGeayan (4.3.9) 
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thesnoustionen beingsdeLined inwSsection 221% 
A numerical example below illustrates the algorithm. 
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Property of Sequence E The transposed inverted E 


sequence is the sequence E being turned upsidedown 


and divided by n. 
% 
(nee Caer (it) i =e (Ares itz) 


Proot: “Srnee from (2.42.6) EEA Ge = E(n). 


Phereroces £ COM 4.30) i, 


es % 
(e71 (ny) "= ce’ (jp! (n)/n) © =I'(n)E(n)/n = E(n)/n . 
(Are 3) 
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x, aw Gl So JEP PAR SNS Sa ie 
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Tiesy(imebere l1seushally referred to as a Fourier 
elanetorm matrix. the sequence Pecan be formally 
defined as follows: 

Definition | [he Fouriter transform sequence Hees 


emoceguenceyGeneteted. by 


Py) Fey ltaty 


Bcnct 1h) = 
af 
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Ey) = sgl 
X, (A) % rwetl ee ee rE 
xe 
x, (n") a rwetl went] | ynintl), ; 


As far as inversion is concerned we can expect 
Ft (n) tovberdexsived from (4.274). But al still faster 
formula follows the definition of an inverse Fourier 
BBddotOrlmmeNOWewem rt nds Chie yinVersesorel (IleLOo be 


-l -2 —n 
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Fo(n) = =| -3 -6 | -3n| = F*(n)/n . (4.4.3) 
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k=1 
(4.4.4) 
7 ; Uk (u-v) no ton ea usV 
k=1 OmerOorea ll mus mcancenwW 2S cy orice, 
Therefore F (n) F* (n) SS Meeps Oe F + (n) = es Gan by fave 
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necessary to conjugate the whole sequence ean tos da — 
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Pomc late Tate Dorie mand ae sequences are symmetri- 
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V. COMPOSITE SEQUENCES 


5.1 The H Sequence 


Definition Sequence H is a sequence formed from 


Sequence —— tas tollows: 





(Him) O (eee) H(n) H(n) | 
H(2n) = = (adie) 
| O EAC ay) er (7) Ett) ei) 
Witch G) en ieander (ede unecde ine (4isoe 2.)a. 
Properties 


1. The H Sequence consists only of matrices of 
the size n = en BeDeGi nga positive: integer. The first 


four matrices in Sequence H are listed below: 
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finn = 4 ee 
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0 


Hneeoremeoece, Lhewproduct of Sequence H and its 


transpose is equal to the Sequence A multiplied by n: 
i 
Hsin) egy yea. Ann) 2 ee C(Sieele.67) 
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Basten (oor cnm sr romes (5. 176) 


ea) Aaa ae 
Buca 2s LO) shows that ay + (n) = A(n), therefore 
ar Lia) = GP CiNee ae (313) 


The inverse of H sequence is a transposed H 
sequence multiplied by Sequence A and divided by n. 

Since multiplying by Sequence A means permuting 
some of the matrix columns of rites inverting Sequence H 
PSmLOentMcaimtOeetCansoosing.H, dividing by nN, and then 
permuting some of the columns. 

Starting from the next section, the operator ® is 
to be employed to simplify our statements for the remain- 
ing of the thesis. A ® B means a Kronecker product 
performed with the two matrices A and B. In Appendix I, 
the properties of a Kronecker operator is described for 


easy reference. 


5.2, Kronecker, Sequence 


There are some sequences which appear difficult 
to invert when taken as a whole, but an analytical study 
of the sequence may disclose that it is composed from 
several simpler sequences. In such a case it may be 
possible to invert the component sequences one by one 


as partial solutions, then piece together the partial 
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solutions in some way to obtain the inverse of the 
OrrgiNal acquence. |The rest of this chapter will give 
examples which illustrate the method. The next sections 
will introduce sequences L and W. These are the 
sequences to appear rather frequently a Pane Two of 
the thesis: The general’ form of the sequences will be 
studied here. 

Let us consider a sequence whose member matrix takes 


ENS Lome scomccction 2.1) for notation Gf Subscripts) 
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i. ie Ae ee ea 
OM OM eee 
Hae) 4a Lee 
“Geib tes til alee) Pete ees 
(6.2.2) 
cele tal abel Seer. Pa hey 
S Ss S sleeve S 
$(J,) = Diy Sor Pee Dy? 


SD Moe Oo ens wee: MO so 


then the sequence described in (5.2.1) can also be 
defined by the 
Detin Perouse ocducice: Ll is aekronecker product 


of two sequences R and S. 


L(JjJ,) = S(J5) 3) R(J,) ; (Seo) 


According to the properties given in Appendix I, 


we have 
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— acl - 
=55 (J) @ R D (See) 


Ease Algorithm The inverse of an L sequence can 
besobtained in two»steps: "first, invert the facto 
MacLices Rvand=S andividually; second, form the Kronecker 
Productsot therinvyerted R and SS. ihe Kronecker procue. 


thus obtained is the inverted L sequence. 





(2.2, 4) Brequires, only igs mubeipiications, plus 


the necessary operations to invert R and S individually. 


5.3 Row-wise Kronecker Sequence 


fetwa sequence W be constructed in the following 
IManner: 


VLU COlSi sts OnlLyeor those mata ces Win) with 
n= JiI5 Ji 15 being positive integers. 


2. The elements of W(n) contain only power pro- 


ducts of two complex sets of numbers t and z, where 
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LOLLOW tie Notations stated in Section 2.1. 





3-eelhe t has onesvaluice in each row of W(J,J,) 
but the value varies from row to row. 

Af. The z has the same value in each row and 
also the same value in consecutive sets of J. rows. 


Uhemvalue Of ez changes arter everye. rows, 


2 
5. ine™=powers are so arranged that the powers 


COietwareukeptsat One fOrsthe =rirst J, columns, and 


1 
increase by one every Jy columns. The powers of z 
increment by one from column to column until the powers 
of z reach a maximum of Jay Phene che spOwerEs FOZ Stake 
from one again and continue to increment by one from 
columm=togeo lLumn, 

Rca ce Meola yeOLMm (2. o 1.) eC leclLOsesmthat tie 
above descriptions can be simplified to the 

Definition Sequence W is a sequence consisting 


of row-wise Kronecker products of sequence Z and se- 
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Yeror Ji55 elements. These elements are split up into Jy 
groups of J5 elements each. The first group only relates 
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[I (J,) @ Z,(J,)] is a checker-board matrix; each row 
(i.e. each Z,(J,)) operates on a different set of ele- 
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Generalization T and Z were defined in (5.3.3) 
and 45.3.4) eas, Vandermonde matrices for4»llustration 
On Vie Reader samaVyemavesnoticedsthatain our proof of 
the fast algorithm, T and Z were treated as general 
Matrices, and that the power series Re ic of 
Vandermonde matrices were not needed throughout our 
jemgiuvabion. elhnerefore the restrictions (5.3.3) and 
(5.3.4) may be removed in future reference to Sequence 
Wo thate@isato say: as far as W is defined as (573.2) 
and T and Z4 are square and regular matrices, the inverse 
OcaWewil bebestiiaterie (5.3.60). 1rrespective of the struc— 
PucCeeO Gel Pandey 

To illustrate this point, we give below 
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Inspection shows that (1) the upper half of W(8) 
has its even columns identical to its odd columns; (2) 
in the ee half of W(8), the even columns double 
their preceding columns; and (3) the upper half of 
W(8) has its odd columns identical to the odd columns 


PHieeiemlowerehal: of Wis). Thus teadssus to rewrite 


wW(8) as 
TCAOe| | o 
W(8) = (Sere) 
mE) Ga 2a) 
and 
he R350 3 9-05 
-.6 Ae .6 wit 
T(4) = - (CASHES) 


We now form Z2(2) according to the two row vectors 


ipa Ase state he 


eT SI | 
Z(2) = : (5438719) 
2 


Z(2) and T(4) are not Vandermonde matrices. Invert- 


ing Z and T by usual methods yields 


1 2 -1 
maa? ae (553720) 
-1 1 


and 
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HYEORR2: OF 
“a 2RETEO: CO 
T ~(4) = : (553 cn) 
oy pik ey oak 
Ojiel =6 6O 
From (5.3.6), Ware e) canape Weomputed sing (5.3.20) and 
Cehs Shs a. 
) 
2 all 
wot(g) = {e7t(4) o| | pt (4) 0| 
-l i 
2 O 4 aah a2 a). 
-1 O -2 -1l i O 2 il 
4 2 O Ose 2 2 1 O O 
-2 -l O O 2 AL O O 
2 25522) 
O 2 6 2 O -l1 <3 -l1 
Gy =) -3 -1 O ik 3 il 
O 21 2) £0 O -1l 6 O 
Oo -1l 6 O O 6 O 
Deel SmeASveLOmVC RUG ya we = if lehie abdieeiyysberer (C5) 52) Ney) 
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Vii wo. ovEMo. OF ROUATIONS WILT 


KRONECKER COEFFICIENT MATRICES 


Ge eeoystem Golution by Matrix Inversion 


Consider the system of equations 
[S (J) 8 R(J,) ]X(Jj5,") = ¥(JjI5,") : (Ones 1) 


The technique developed in (5.2.4) shows that (6.1.1) 


has the solution 
hi Toys Wer Ge ee Rac wen (Ee) 
1 2 ik leet oa 


al : 
(35) 


Suppose R and S belong to the EF sequence so that R- 
and Smel5) are simply Ria on and S(55)/J5 (see 


Ae Oe ee en) 
X(J,J5")= [S(J,) @ R(J,)1¥(5,55')/ (545). Cha eee) 


Now we estimate the number of operations needed 
EO GecCOMmpUcd NO Nr. 


TOMS GLOLMetiesaArLLHMeL2CMOlmn Onl. 5) —requaiies 


oe D) ei: - r : 
Ji55 iisiet ec Once Omyue Us tile Jj55%*J145 matrix 
be Gh 
(yd) e More multiplications and additions) to ger 
[S @ Rly, 


JiJ5 divesions to obtain’ x, 
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Additions usually need much less time than multipli- 
cations. The CPU time required for additions may be 


ignored for the sake of simplicity. 


. oe 2 oy 
Total operations= 2(5,J55) ap J555 ~ 2(5,J5) (Get e4 ) 


oat Ji 15 are large. 


This indicates that using matrix inversion to solve 
(6.1.1) has a distinct advantage over using Gaussian 
Elimination, which normally needs ney operations. 
However, (6.1.3) does not give the quickest solution 
HOmn(O eb) melee a Swamstiull sbastersalgori thm based 
on transformation of the system (see 6.2.la). This 
algortunmecanecolive (6e40a1) fin J,55(5,+5,) 


operations. 


6. Zee COM LOnmOt eaa2-DlexecvSctem byekectangularizing 


the Unknown Vector 


Definition A vector X(n') is said to be rectan- 
Gulanased ito aemMatiax A (/h, hk) ots. NeeleMmentssare 
SOuUiminCOmkeOLOUDS mL athe sOniGidld morder, and mie 
groups are displaced sidewise in k columns with n/k 
elements in each column a. An operator Rect, is used 


k 


for this purpose.. For example, 





* 
Note: We assume n divisible by k. 
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If X is a matrix, then to rectangularize X would mean 


to rectangularize each of the columns of X. 


x) X, X72 Xj9 ++ eA 
Rect. (Rect ,X) = Xo Xe Xo X44 ree X53 
X_ X¢ Xqg Xj0 «> <5,| 

In general, Rect, X (a,b) =e (a/ky DK. 


Definitien’ A linear system of equations Li(n)xX(n')= 
Me(Nee) sus) Sardmto be m-plex Wi iits coefrLicient matrix L(n) 
is a Kronecker product of m matrices: L(n) = S, (Jy) 8 
S.(J5) Orr atete AO ee UE, n= JjJ553---d,- 

Theorem 6.2 Provided L(J.,.J.) = S (J) @ R(Jj), 

joes 4 + f : i] = q 
the 2-plex system of equations LiJjJ5)X(Jj,J5) Y(JjJ,) 


can be transformed into the following: 


R(J,) (Rechts. x)= (Rect ay) foes jis (6221) 


wherel(sec eSectlon 2ell for notation sos subscripts) 





Note: Ekstrom (1973) defined a "stacking operator" to 


perform the reverse of a rectangularization process. 


BA 







: - t s * ‘ - 
clung oiet o4 dethh «treat ih areas = 
; E 4 } BS 88 eee | oF : - 


> aera a ee eee = MS + fay 
%% t * if “ert . © ; ig Tae ats : 
is t | 5 


Ooo. ee itn ey - (. Che ine Bsn ivs<) t') 


‘-. r ae | a : ¥ ,. het po 75? Fi sits 7 n nN : : / 
orn? . i has Oh Taupe 4 teint) —= 
7 ott 4 Lig? 343 aon | i, ob Remaee ae : 


1p pee dine iat. 1 youn 
(oh An 2<% Re atk Ae 


are oA 





) 
Sie! Soe oe SO: ee oe altel 


xX, x Soom os x 
a 2 ple te (I2—=1) J 12 (6.2.2) 
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and Rect_,Y is same as Rect..,X except that x is replaced 


by y. 


J2 J2 


In ouner words, 66.2.1) dives theysolution of the 


2-plex system as (6.2.1a) 


Sa 


= pa 


~ agen 


Soir 


pe 


Rect_.X = [Rot 


acl 
~ Vi Stesry ie Mae #2 la) 


(J,)] (Rect, 
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(67.253) 


Puteinge uc Ist, (J, +1) th, (25, +1) th, (35, +1) th 


equations together and combining r with x to form a 


term, 


(6.2.3) becomes 
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= | (6.2.5) 
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Then, we put the 2nd, (J,+2) th, (23,+2)th, eee eee 
SQquacrons sO tn (Ones) Ss COdeCtIer,, endsinethie Same manner 


we shall obtain 


ve aah es 


‘ait 
(Guna) 


Ie ie 


on sh ios 
Ro (Fq) (Xsy yy e+ Xozq) = SQ (Sa) (Yo V5q40°°*¥(g2-1) 5142) ° 


CSAs US) 


We can also put the 3rd, (J, +3) th, etc. equations 
Ofs(G.2.3) togethereandrobtain another set or equations. 

THeBneX<testepewtlie bewtotputm (G.226a) 7 (6 S2Ki/a) 
and the first equations of all other sets of equations 


together. This will yield 
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LONG: -l,—,,t 
= Rect,,Y[S, (J,)] 3 (On2 20) 
Prem@ntcap lyingetoe2eo  ewieh Rew (as), 
oo eel = eet 
(x). 22X57) = 18 (J,) Rect,,Y[S, (J5)] 4 (G2. 9) 


Seite Given Uinta 6.2 0) ,mel Oe 22) sere. Logerner 


will give 


1 j=, ,t 


-p tl vIs7 
=R “(J,) Rect ,,Y{S,°(J,)]. (6.2.10) 


(x iy aoe a 
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Comp nG mG ) pe Ol LU) andasimitar formudas ,. wild 
give (6.2. 18). 

Fast Algorithm Theorem 6.2 shows a way to trans- 
form one linear system into another as stated in (6.2.l1la). 
Lie tilesCOULSC OO metic one co rmien the column vector 
variables X and Y are rectangularized, At the same 
time, the coefficient matrix L is broken up into two 
small matrices: the premultiplying matrix a and the 


its 
i . A necessary and suffi- 


postmultiplying matrix (S™ 
Clent CcOndLUions TOrithe transtormationeiS Us =15 GO 5R. 

The speed of computation can be further improved 
ii R and S belong to the Sequence BE. In this case, 


-—1 “a al. 1 (NY) 
R a R/J 4 ands) (ole) a= S/J5- Therefore 
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REeCto wes R(J,) Rect 12 


J2 ce, 


Then solution of the 2-plex system turns out to be a 
MULerpuiCattoneorethiree Matrices, and aldivision byin. 
The total number of arithmetic operations is 


2 2 = 
Le Ue a) acre! J5 ~ Jj 55 (Jj, td 


12 199 7 sie. gel J>5 are large. 


2 1! 
(Ope 12) 
As we all know, Gauss Elimination is the most 
popular method for solving a system of equations with 
non-diagonally dominant coefficients. Gauss Elimination 


requires ise Munem On muta plicacvons ian this case. 


i: 


Comparison between these methods is illustrated in the 


following table: 







Estimated CPU Time on 
IBM 360/67 machine 


(b) 
Fast Algor. 


Number of Complex 


Poe 


(a) / (db) 


(a) 


Goussmo ain. 


Unknowns n=J,*J5 















55 m-sec. 5 m-sec. 





100 LASSE Cre 80 m-sec. 166 







L625 min. Le SeC. 
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6.3 Solution of an m-plex System by Rectangularizing 


the Unknown Vector 


The technique developed in the last section will 
be extended to an m-plex system L(n)X(n') = Y(n'), where 


inf) iSeaskKroneckers product. oL am Matxices: 


ITs yee Jege—ese (37) 8 Seo (oe 5) °e..-8 So (o,) S¥si (5). 
(Ono ol) 
THaecOocemeaGeoe elt tn) = S (JL) @ See Wa) a 
system L X = Y can be transformed into 
oat ety ae 
S(J,) Rect, ,X = Rect, ¥ S--1 e-1). ° (One) 


PYrootsseLromes (Oe ee) sutease clear that) the solution 
OF 
Sy ® S x = Ve (Goa) 


Ls 


1 mel ass 
Rect, X = ney Rect, ¥ S. * (Gnoo 4) 


Transposing both sidestol (6.3.4) 








. P< 4 ¢ “ — ~~ 
nis ; J ; aS BiA @ pid Cedi 8.0 Se 
et - — ae Gat & 7 tied ay tae 


meteve setken ae od See al 


4 ern -4 ales al ih Lj oie hes at oe = @) fata > ae : , 


j 
7% 
- fr 
[. 
(oy a) fe i. 
_ nt 
iy ti £ f Jess ILS 


SS ee 
A hoee 


Cn oe : 
_ 
iy ye ca 
: i. 
7 
- ea 


+e > Fi 


(E-2.a) aw 





79 


Ce eo. -t,t .-1 re - 
eee ei ee ee ome oe a Rect OY oF 


ites 

-1° (6% 3,5) 

Premultiplying (6.3.5) with Seecotabddches | (670.2)— 
Theonenmows ames the Coectii cient. matrix of an 


m-plex system takes the form of (6.3.1), then the system 


Cans be geransrormmeqn) 1nto 


t e t t A 
Sma STE Geel) EOS ie hy 0 OSI Te Rect,, X = 
t : t pa Ee 
Reet yma 1) RoC ts (moa) $22: (Reet G84 )S, Jiveee: vee) 
—t -t 
ie) at sc (Ome) 
Proof: Let the system be 
S © S- © gsq By Oy =e | (Geese) 
TnheoOcemab sol ewil le transtor §(6.3940)) ante 
Ss @Ss @ a) uy ae Shae we ee (6.3.9) 
a aay Sake > Rect,,X = Rects, Ll: Bae 
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Rect 5, (Rect, YS, )S, 4 Cae She, 85 OP) 


The process Gontinued in this manner will’ y¥eld (6.3.7). 
(6.3.7) VYeads* to the 
Fast Algorithm The unique solution to the m-plex 


system 
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SO Sh aad Sey oe es (omeo 130) 
Ls 
RECES (yy) RECES (mu7) RECES (a3) 1 *RG2REC Ty aS 
she ecie ey (Rect ma, (ns (Rect, (Rect, Ms 
Se cae Se JS oun. Cane 


6.4 Matrix Inversion by Rectangularization Techniques 


Finding a matrix inverse by solving a system of 
equations is a well-known method. Now that we have a 
fast method to solve a system of equations, we hope that 
our fast method can help to invert a matrix. 

heteuceneqimanyereplacingi the data lw int. 2Z.)) 


Wa Cima sUll Ivan Me tick ox) ky (I) Ps 
Milt ab (tee Ln) (6.4.24) 


Solving this system of equations, we shall obtain X(n) 


WiLehmeLS@Meniecmtniverse Oem (nN): 


Se Saray ae Saree (6.4.2) 


In a special case when M(n) = L(n) = S(J,) ) R(Jj), 


n = JiJor ehe=system (654.1) can be solved=by thesrast 


method (6.2.la), namely 


ib i 


Rect (nh) =e Ra (J,) (Rect cas (6.4.3) 


g2'h 72 


P 
s¢ te / 
a 
bel 
y NS ' 
“eo 


(tr had) | a (Age Sees ® 















-< Mi g'din ch nie i. 


i= ai 


d - ay - 





‘ ; 
y 4 D lay te 


“ , Z / 4 , i _- 
=} ae “ 93 
oy) ; i | j 4 mn ve ae ta» ’ {i ‘see 


Lupipamed. Yel ase Siu ag 1p 
t tS a's fs sniiwiot : Mon romeshig LP ie Beale of pagans 
ie? Peeaee Fait, Ren art E SW. ad anti Final re 


i Pr 
ere 
, 
sqpAa ow “tos i407 Tee Meee S184 oe i See eae Jest : 
Tae (sal ue aber oti es. wae 
i's hes Y, 28 +12 A riselqas MG <i aK? nk oe _ ; 


iisbt xeon oe ae a, ¥4 
what) (yt 23 aidan’ OR 
r 186) Lost ae mesa coaster 


> ae adrqval ast ed a bei 


1 pyle ic 
(S.. 4.9) | woes = et es 





% iene = (she - (sty 


aoe i ; : 
4 = , r on h 
ar | 


- 
7 an aa 
ier 7 





2 


Then we have the in lea after rearranging the elements 
of Rete, (Le tl) Since Rect,.1 Contains nothing but 
zeros and ones, it needs no arithmetic operations when 
multiplied by another matrix. The whole equation 

(or4s3) needs as MUubCEpIacatvons. “This 1s “the -same 


number Of operations required for - calculating eee 


Wienkate? WSR eet) un 


6.5 Summary 


We provide a table to compare ail the methods 
introduced in this chapter. This table does not include 


the necessary number of, operations to invert Rand S. 
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AIM MATRIX INVERSION L_ 


Kronecker 
Gauss elim naw Mat. onws: 


5 
METHOD PasemeALgoi. 


FORMULA 






a 


ecuine = 


Rol (Rect 16) rc 


NO. OF 


MULTIP. aes 
REQ'D. 
REF. NSS GES (Gia) 


(1968) 





AIM SOLVING A SYSTEM OF EQUATIONS LX=(S®R) X=Y 





: Kronecker 
METHOD i | Veyos aeoee Bast Algor. 


i 


X=tor ¥ X= (s-teR7 


FORMULA YY | Rect X=Ri (Rect Y) aoe 








(Ome) ee (Oma sa) 





Westlake 
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In the light of the above, we conclude that 
(Ome oy ei smenem oes tealgori chm Loresol Vang a System 


OnecdUatronss With le— ro OO Reas its cochiicient matrix. 


However, if our aim is merely to find an inverse of L, 


Eneimoees manda (o.4.5)sare equallyseiiicient, 
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Vide sOURTAPPROACH TO A COMPUTATIONAL SOLUTION 


(lee hive: Physical, Problem 
PatoameCCalaAiwl reatment ofea Vector Field 


A time-harmonic of an electromagnetic field has 
two components: an electric field and a magnetic field, 
each with three spatial components. Thus the electric 


Fieldsis welttensas: & = (ih oF Hane and the magnetic 


field as H = (HLH H ee (iE SmOULBUSUaLeDLactiCe: LO 


Va 2 
Studveonl yates directional component of thevelectric 
field, which we call Eo. The other components will 
behave in a similar manner as EL 
The component EY GomomcunCctTONnwot times and space, 
with a frequency w and a leading phase angle 9. 
Mathematically, it is 


EL = E,(X,;¥,Z,t) = Pa Vir 2) COS) Eto Gupy 7 2) 


SRE (U(K, V7, 2) CxD\— LUE), 
where u(x,y,z) = A(x,y,z)exp(-i¢) is the complex ampii- 


EUCemMO Leis. 
x 


TMomSLUCy a hemopt iCal eficld ites cUaerc Pol icG 


Sctudcyausonlys 


7.LbD, Propagation in tsotropic Non=homogeneous 


Medium 
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An optical wave in vacuum satisfies the homo- 
geneous Helmholtz equation. In the presence of a 
semi-transparent object with refractive index n(x,y,2Z), 


which varies slowly in space such that | (2n7) /n? | << (10) 


14 
meter - (seemAppendix Vj, the optical wave will be 


propagated according to 


Foe c= wae. Gigi 
Equation) (7. lel) isa Sehrodinger's equation, where 
Ko = 21/A is the wave number of a coherent monochromatic 


SOUrCEE Ore light, «and va is the Laplace operator. 


(elcome liemoCateered sa2eld 


Ine thesspaccmwiere Light Ssescatlercd, tieeoptrcal 
fHeldm@rsecomposea Oretwostsaelds: the Ssource-free back=- 
GrOundmiret| dmancmrnecerteld Scattered by tnexob ject. 


hetrwussaenote these fields with Us and ME a 


(e) Fy a ; (jee) 


Since Us propagates in a homogeneous medium, it 


can be described by a Helmholtz equation: 


Wo ie a Ste re IE ee a ake 
ae fey yoke) ! Oo 6) Sey at a 


SUDStmUULINGMN (Ll. cyimanden /el.s)into (fol. 1) yields 
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if 


2 2 2 2 
+ =- = 
V un KS es KO (neal) sus. (ToL. 4) 
Let us define the scattering potential ae = -k* (n*- ik 
Then 
2 2 s 
V de + KS ‘eae f he Cia sy) 


Thais (eleven Sy eoeche Approximation 


im Many Cases ,, the field aes scattered by a semi- 
transparent object is weak compared to the background 
field Uo: Rome Mit yen (el >) wemsiall, saCCording co 
Born's theory, replace u with us, so [elgvelias | Telele) va Wolenel cgysheve: 
SrdemO tees becomes Independent, OF aoe Then the 


equation for the propagation of tte 


Vee UL +k” u =f U (ciecwlen tO.) 


is a non-homogeneous Helmholtz partial differential 
equation. 
Pepara CuiacesolutiOleOl Meh. .0 sl Sm (cee Appendix 
VI) 
u = | Geo Ue CV. (Le) 
sc siteo 
V 
with V representing the volume of the scattering object, 


aV =dx dy dz (see Fig.l), 
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G denoting the Green's function given by 


- exp (i 1 | = xaly 


a a Ve Sah Gieatste 


—— [x,y,z] ° being a space vector pointing to the 
DOsTeroOneOteGVeInsidestiesscactering object, 

X'= Ix',y',2'}° being another vector pointing to 
ae DOUCt Cow i, Ze) Cie tei etieostar—t el 
region where Lee is measured, or within the 


scatterer, 


7.le The Far-Field Approximation 


Subject to the condition (see Schmidt-Weinmar, 


Lo 7) 
2 ' 
ROR. <span e a )) 
| “| tan | XeX a 
|x'- X| can be approximated by [X'| - Txi7 for the 
numerator an s(7.1.8)) and by ||X"| for the denominator. 


TiemEOrmcimweapproximeatton ac 1llustrated an Figure 2s 
Heres" -! indicates an inner vector product. 

We shall introduce a new vector K whose magnitude 
is always a constant equal to Kor and whose direction 


Will bewthe same-as the £Lar-field pointer xX!: 


(K Keen Kg) ae (75110) 


K =k ’ 
x ys 
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exp (ik |X" |) 


us = | SEO GOt SONY) (eee )eatehien is ea) 
sic an|x"| Vy sO 


ToOwCOnvert the above to a digitized form, we let 


Bi e8 4 EL 
Fo F (Cs ees) 
exp (ik, |X'|) 9 (K) AV 


where @(K) is a form factor, and AV = AxAyAz is the 


digital equivalent of dV. We also let 
fe — fou 3 (CUA abe) 


Appendix IV gives details of the mathematical derivation. 


A result is yielded as 
Ea et) ee meer: (he)! fac. CES) 
Me Gye Fe 


DLeesioulisbesnOteOatiat slearepresents sthe ssounce 
densities due to the multiple Se ng sources 
within a volume element AV that has its center at x. 

Lee le ee Ser DCOCUIC ta Of uy and i both being 
functions of x. Hence £ depends on x. f does not 
depend on K. 

KetsSeanveclLoLeoucOns Udllem Leng ul Kor and aligned 
Timechemd! rection On Xe (see e/a 0) te spO1n Gu 
moves about in the far-field region, K sweeps a spheri- 


Gal Surface. Any point on this sphere, as represented 
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Dye pehasmoniy @twa degrees of freedom instead of 


three. p Lf WS and KY are the two free dimensions, 
then the third dimension is determined by S and Ke 
Hence K = (K_,K_,K »E: K_ = Ke - Ke = ae 

Sai 2 z O Be y 


| ae meeport ronal to ae (SESe Ppa kas ihe 
represents the iscattered field in the far-field 
region. The magnitude and phase of Us, can be 
measured under far-field conditions (Schmidt-Weinmar, 
1975(1)) and hence F may now be considered as a set 
OPaknownuequantities, which are the input data to our 
inverse scattering problem. Since ug. Sma IC Lacon 
Cres awirchievcmamhunction OLR, we concludesthal F 
has two degrees of freedom represented by Ky aueyek 1 F 
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7.2 Discrete Representations 


Tiemrest@eonkcausmcnapler WwiElilgbe devoted to 
defining our problem in terms amenable to numerical 
analysis. This section describes a discrete coor- 
dinate system in the space domain. . The next section 
will establish a system of equations suitable for a 
computer solution. 

At the moment, we are dealing with a semi- 
transparent weakly scattering object of a finite size. 
We shall assume that itis contained in a rectangular 


block-shaped space. The space is subdivided into 
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rectangular blocklets of equal sizes. The scattered 

freld source wretnin each blocklet wil be “averaged 

endeacsuMecmeOnoc Ormonunated irom tne blockilet center. 
Our spatial coordinates are chosen such that 

the entire scatterer ives ne ties ti mei soctamt,. and 

axes are parallel to the edges of the block. If we 

call thekdimenstons of veach blocklet Ax, Ay, Az, aes 

the coordinate origin should be placed at a distance 

Pax Ay, Az]© Frome the centereor a comicreblockiler 

UE eis Sh 
jivets, cl) J 


Ue be che numMberneon blockletse im the 


UL EG ae 


X= 7 ee -OEeCUNONS.6 LGt) also Iye Jor 3 be the count 


numbers in the three directions. Let again X=[x y z}° 


bewthe Coordi natesPoLtstnercentersmotecie blocklets. 


Then our choice of coordinate system above results in 


x = 5, 4x Jpet ls 2p eee rd yr eee TZ) 
y = j,Ay Jn={1,2,---rdogyre ee ITa} (F200) 
z= j,Az2 J3={1,2,---r3357-++ 153} 


For convenience in writing, we shall label the 


blocklets with an integer v where 
VE Fy Iz Goyrh) +5452 53,71) - Wis. 2) 


(7.2.2) shows a one-to-one mapping between the count 


numbers aioe aig yew al ee and the new label v. Obviously, 
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A SCATTERER AND THE SPACE IS SUBDIVIDED 


ifs? OCCURIED INTO BLOCKLETS 
BLOCK-SHAPED SPACE 
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THE vth BLOCKLET 
IS REPRESENTED 
BY A SOURCE POINT 
fy AT ITS CENTER 


(Xy Yy, Zy) 





jo ae ER ORE 
THE SCATTERER IS REPRESENTED BY J, Jou 3 SOURCE POINTS 
(The numbers attached to source points are the values of v) 


FIG. 3 
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eon mans, e202) Vos {1,2,...,5,55J3}. 
RL ievagpiablLes assoclated with a bDlocklet can 
Mowebes label med lhe first Varlvable involved is xX. 


Pica chiim Gd svacouexeeanademodiiyings (7.2.2) , 
KOS [x zo = (4, Ax, 3 Peis ene) 
vo vilyt Vv = laa! ' JoyAVr Ja5, 42 e ofe07 


The second variable associated with blocklets is 
feo S eae pot ner! Calwcuenclty,. Located, at the: center 
of each blocklet (see Fig. 3), to represent a mean value 
of scattered-field source density within the blocklet. 
TiewOctin bev onmote: slemg1Venuby scalaci1on (7.1.13) as 
bbe producto te. the scattering potential and the back- 
GEOundecveld ime S ea LunctLOn OL xX. We shall attach 


ehewsame subscript v to xX and £: 
f= tL when era keel Ci.2..4) 


Since all the £ are subscripted, these te LOD — 
1,2,6++75,5553, will be listed together to form a 


VECLomel (it wawlLel dg) seu 


ed Pest dee See ee is (7.2.5) 


ee 2 gig2g3! 


Now the space domain has been discretized to form 


MUunELOrmily spaced Gria Structure.. "Vhegurad points are 
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finite in number, labelled from 1 to J,J553, 


located at x, where x Pomc inied witian ey) = 2,.o eae The 


and 


scatterer is represented by the grid structure in 
WorchmlLoreachneqrid point we “attach a hypothetical 
QUentity called the scattered-fLield source density 


ty WinCE SeOein ined sii (6/7 2. oS) a. 


ie eelner Sysrem of BQuations 


Bach spoint on the grid structure, except those 
points with ft. MOP SmCOnSTOCreCOmt em LOCetLon of “an 
ince pendentesousces Of Optical disturbance. It emits 
an optical wave which propagates in a manner described 
DyecdUaeloum(/ se lone (hescombined) eltcet of all sthese 
SOULCeSMCLreatesmtnes actual optical) field sin the far 
fleld region. 

Wiesscabterco eral ftield aserepresented by F. 


According, Copeduacion a(/#l.14)),) Peis related to ft. by 


a ) exp(—ik =x) oe C7eaee) 


In our inverse: scattering probiem, we walit Lo 
determine the unknowns fie for Heke dy £5345233 from the 


data F. It is necessary to select J,J5J3 independent 
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equations to solve for these unknowns simultaneously. 
How to select these equations to expedite our numeri- 
GalmcomputLatton» 1s a topic to be dealt with in 
Chapel ment eis stage Tiemay bersurficient to 
eopnnte that we can choose, in some Suitable way, J,J.,J 


aS) 


vectors K, labelled as Na with 
u= 1,2,-++55,5953 : (URES ss) 


The sets of data F should correspond to Rue ands pe 


labelled as he 
F = F when Vere ee COS SEC9) 
Aerang inGm Lic a imethe order Of U,pwerllave Cheyvector 


BaGrals) 7. 


Rs t 
F(n') = [F, Fo --- Foigog3)° UAE eS) 


Then the selection of io yields a system of 


equations 
aS eee GX) ee 
Vv 
Py = Wexoiaiks x=) eet CUS Ae) 
V 
aralyparcy L SO 0 ee 


Writing (7.3.6) in matrix form with the help Of i. 3. 0) 


Sais) AG Ae, 
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Eto) eae Pn) fren) (ie) 
where 
exp (~ik, +X) exp (1K, -X,) Sac exp (-iK, X54 59353) 
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Piya exp (-1K-X_) - (i322) 


Since P is not diagonally dominant, solution of 
(ino 26 Eby sterative methods wil linet converge. An 
attempt was made to solve (7.3.6) by Gaussian Elimina- 


tion; the computing time taken on IBM/360 was: 


No. of Complex Unknowns Required CPU time in seconds 
25/ 3 
64 | 37 
400 532 
2000 86400 


For large number of complex unknowns, solving 


(7.3.7) by Gaussian Elimination should be impracticable. 
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VIII. METHOD OF SOLUTION 


oe POUELINe=oL the Method 


Oureme rhodes nasolving PE = Fis guided: by the 
following: 


K ee (eee Oy 


a ema We se Ole (he = ar Kae ee, 
u HEA VARI SrA 


ei pila. lene st aken tompe ml dirrerent 
points on a sphere in the K-domain [see 

equation a0) i) en toi tousphem.caly surtace 
is continuous. We can choose any n points 
on this spherical surface as Kj, 1Kor---+sK,- 

2. The selection of KL will not affect the 
unknowns £., because fy Moma uUncoLOn Or Xx 
NO teow Lea (see Section 7.1). 

3. The selection of Ky Wa biabnect ie and as 
(see 7.3.9 and 7.3.2). We hope that a clever 
selection of on MAVeslNp ist yar LO. One Oo Le Oui 
Matrices in Part I. If this can be achieved 
then we shall apply the fast algorithms that 
we have introduced in Chapter II through 
Chapter Vl 

As we noticed in (7.3.9), every element of P 


involves an inner vector product Nee which is a 


sum of three terms. 
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PPevOUlmbemOMEr Cult stor wdenticy Powith the patterns 
of any simple sequence that we previously introduced. 
The best we can hope for is that by a selection of Rue 
P can be turned, into a product of some sort involving 
three factor matrices, so that the syotem PEG=—R wi 
be solved eventually by matrix decomposition techniques. 
Our method is based on treating P as a row-wise 
Kronecker product sequence W, consisting of a 
Kronecker sequence L and a Vandermonde sequence C. 
The matrices in sequence C can readily be inverted as 
a partial solution. Then the rest of the system will 
be handled as a 2-plex system with L = R ® S as its 
coefficient matrix, and R and S will be inverted as 
another twoepaktial solutions. All these partial solu- 
tions are then pieced together to yield a final answer. 
The following sections give a detailed explana- 


VOI. 


8.2 Creation of Row Redundancies’ 


Lethusetrewn ce (7.3.0) inva tormm more Convenient 
LOLeManlpulacwon. aby a GO. .l) aNOee (7.26.5); ) WCU tilt AY 
indicating the integer “part of a real number A, and 
(B) mod(C) indicating the remainder integer resulting 


ELrOMEe, C. 
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elements will have the same factor one 
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Soo0) Creation ot Column Redundancies 


Eng ondem cOelurna (nn) intoWiny, only, redundancy 
in the rows of P(n) is not enough. We must create 
column-wise redundancy (see Section 5.3). 

| Column-wise redundancy can be obtained by select- 
ing Ky in special ways, and there are two possible 
patterns of redundancies. These will be discussed in 
detain Chaprewes ls llmetivs chapter, (6b US assume 
that we can succeed in splitting the n variables ré seh exe) 
JiJ5 groups, and assign one common value to the J3 


‘ q . ‘ 
variables r, in Gach Group, This meens 


ry a req = constant 
Cores sey) 
ay {Int[ (u-1)/J,] }mod (J, ) 7 aw 2 hy actos « 
In a similar manner, we Shall assume that we can 
achieve 


Su = Si = constant 
GOs. 2) 
ky = Int{u-1)/3,5,}+ Dy 2 oo 
(CSeiseal mec cae ors ec) Chang CmGos. 2c) aero 
Syleal? ele fees (aeeee) 


Riva uesklacke su 
We did not assign common values to the variable 


t,- Nor didsweschoandesthe wUbscrapt uu for ti. | fhas is 





i¢ 












- 


7 ey , 7 
F - : 
corti Wie row a0 mE { os 
YVoreaiyoen W Poe = . 
a¢ 
’ wipe Ray 1% 
) AC apR eae sobaithoer: 2 


~~ ; - *5 
~jaolee vd hentsido ad ike Varies oxiesiadie 
afer) i>o) \ 75 “arid Brlh poten ‘d2tbegw oi _ 
at. SSgau9 sch ea tt iw Saat? > ts 2h guy ta nee, 
ater # Biase 


a 


y ry ‘ oat | 
Sat -t 2alee)3 iy BS eae sed ee is pret te ?, 


t>) 


. : i] a 
; j 4 = 7 
ifs 494 Satey (One Tanto? Ae Sek Sak { re qutep, chy - 
‘ » * a } x aay “a : a 
row 2ldt .Auedp/ ies 2) 7 eel 





- 


- Z iy 
SmusSes. 24 jel ,Yeuqgsaag ainda: al 


~ Cae | 
— e 


J 


: : 
LDF ORY 3 . de Gt 4 7 
< . ; ie 7 - 


, rs o 4 in | - _ 
ee eR ae ie 


‘ip tegal = is = ot 


fal ow 5 (se hi f= aw. TeceRem ipl uke 4 ne 








D2 


due to the fact that only two degrees of freedom Ky 
and KY were assigned to K. Ne will be dependent on 
Ky and Sy (see Section 7.1). Therefore ty cannot be 
assigned a common value but depends on the values 
chosen for re and Su 
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Letting u = 1, 273704475, 5553, the whole matrix P can be 
WoleEen OU. 
Subs ec tetomene. cormat (8.0-4)saP (mn) can be re- 


written as 
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8. 5.aComputing, Time 


The approximate number of complex multiplications 


required for our fast algorithm is listed below: 


Main Actions NOwen Oo GEC oll lex 


ee | ipli 
rocedure Required Multiplications 


i ibeherbliee qelekays | Read Re Op, lies een Nil (see Section 


ri scl. say) 
2 LN VeELeL JiJ> Parker's method "Hales ie) Wh 
Vandermonde cae tae ele 
matrices Sate SEEeS ee Gee 
me (53) + 
Be Gompute se! \eMulLtiply Baise) Oh a 
7 aie 
with me (J3) mene 
m= 1peee 15455 
: 2 2 5 
4. Compute £ Mubtcipiyechree J5J5t3 3 I5t3 {5953 = 
Matrices of sizes ; 
J,%T41I4%XIT51 (JT) (J,4+5,+53) 
J5xJ 5 and divide 
result by n. 
The above sum up to I,J, (4.50540, +3543 3) complex 


MU WoLpilacations: 


Let oe Jo 


equations to solve. 


=ua— 20so then wethaveba, systempor 8000 
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Using our fast algorithm, about 744000 multipli- 
cations are needed. If IBM 360/67 machines are used 


for computation, solution of the system of equations 
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needs only a minute execution time. (Compiling time 
is not included in the figure as the program and sub- 
routines can be compiled beforehand and the object 


program is stored). 


GPOmPotlOrage Space 


NoustoOrer a, COCEILICleEnt Matrix P(Jj5553) would 
require GC, DTN complex words, whereas storing its 
Kronecker eractormalLnlicesen, oO, .anG me only needs 
Tends 4055505 COnplecawoLds me ASmemiict@ici OlLeLact, 
we need not store the whole Vandermonde matrix. We 
can store only the first columns of R, S, T; then the 
other columns can be generated from the elements of 
Chelmer rote Co lM ne sO sdOunguwemcan reduce tne 
storage space to Jit Jot J55553 complex words only. 

In the latter case, to generate the other columns implies 
extra multiplications. Whether this is worthwhile 
depends on the size of T(J HES J is big, the 
lattertis preferable. 

Let us assume Ji= J5= J3= 205 MSeOrI NG atilemmat ta x 
P normally requires 64m complex words. Storing R, S, 

T in full needs 160K complex words only. Our algorithm 
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generate the full matrices R, S, and T. This extra 
HuMnbereerMecitceecase 15) 30000 multiplications, about 
10% of the 744,000 multipild cations gorigunally required 
DuweOULm Bact wal or acelin. 

Furthermore, our algorithm handles the Ji55 
Mates ces me! one at a time. A large number of work 
Spaces can be assigned to the same core space. Thus 


it reduces the amount of core storage considerably. 





Dx oh LeCl LOND OF —PHASE ANGLES 


We showed in Chapter VIII that a fast algorithm 
can be developed if we make some of the sale take a 


common value r ang also make some of the 36 take a 


kl 
common value Sio° We have, however, left out 2 impor- 
tant points unexplained in Section 8.3: How can we 
make ry and Sy take common values? What is the phy- 
sical meaning of equalizing some of these parameters? 


The current chapter will answer these questions in 


detail. 


OMigeeP rnc lpLesmoOneoeLeCctlon 


Theselements of P are given by (7.3.9): 


Davie exp (-ik,-X)) . (9 1a) 


Kk, may be considered as a phase angle of Puy the 
phase angle is resolved into 3 phase angle components 
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The relations between the Dee components and phase 


angle components are stated in (8.2.2): 


sets exp (-1AxK, |) 

Uy nll z , 
ee exp ( iAyK,,) : (OR a) 
ter exp (-iAzK,) 


In order to have repeated values for r! and s!, 
Woicheares per odlesrLunce1ons, 12S Only necessary to 
select values modulo 27 of AxK and AYR OrweCu eho 
rent values of u. 

NowaweunaVemGUberlcs Guldingeprineciple for 
selecting the phase angles: 

We must create redundancies by selecting values 
{NGelsiey Aap Gaeore AXK Ly and AYK Vu: 

We note that selecting a common value for an may 
WesuLeminmae sing larematcix eh (Scemo.4. 9) ewhich contains 
only rp, as its elements. Similar problems may also 
arise when we equalize the St values. Furthermore, 
the selection of AXK, and AyK ou wild automatically. 
produce the AZK values (see Section 9.2). AZK, Wace 
vig" turn, determine the values of Cee and hence the ma- 
trices me! which may bewsingular asuwells  “Theretore 
CUEESeCOndepriInclples1s: 

The selection of phase components AXK and AyK ou 
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coneerned wills be. infeasible and our fast. algorithm 
Witleleeiatel 

Besides these principles, we have our third 
Principle; namely; to select the phase angles so as 
to stabilize the Vandermonde matrices me* Our 
fourth principle is to select phase angles to turn 
matrices R and S into E sequences, etc. All these 
principles will! be explained in the subsequent 


sections of this chapter. 


9.2 Graphical Representation of Phase Angles 


Heteicgbescemembercdythat Keeista three-dimen- 
Sional vector meta eee two independent components as 
it is chosen on the surface of a sphere (Section 7.1 
refers). If we select ee and Kou LOmmonDaknt citar 
u, then on is determined by the equation of the sphere. 
Since six, Ay, eAzeareeGhosen bysoursphysicalmprobiemegand 
we are not free to assign a value just BOreENewesakero er 
derastealggmalhim mawe Canmonlyeselect. twoseolbuthestiree 
phase components of Diners Usually wesseitece AxK, and 
AyK and leave AzK oy to be determined later by the 


selected values of AxK,, AyK anda hz. 
% yu 
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In other words, we can now simplify our selection 
of phase angles to the selection of points in the plane 
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required to select n points (AXK, ar AyK,,) nike 4st ale ls 


plane. These points are actually the projections of 


n 3-dimensional phase angle vectors (AXKL AYK as 


ea) onto the (AxK se plane. 
ine Aas and AyK aLessubyecemeO nO srestric— 
xu yu 
tions, the n points or phase angles may be randomly 
scattered on the (AxK,, ES plane as shown in Fig. 4. 


In a real-imaginary graphical representation of complex 


numbers, a and st Esrminacte Ol a Unit circle. tare and 


s) take phase angles of arbitrary magnitude as ailustrated 


iM Chem bOtLOMeptCLuLe Of Fig. 4. 
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value Sp5 (see 8.3.2). 
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rows identical, and consequently becomes 
singular (see 8.4.3 for ,T(J3)). This 


requirement means 
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(2) Interchange of absolute values of KS and K 
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two 45° inclined lines. Then the vertices are the 
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FIG.6 THE FORBIDDEN OCTAGON VERTICES ON (K, K,) PLANE 





FIG.7 THE FORBIDDEN ELLIPSES ON (AxK, AyK,) PLANE 
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forbidden points for oe in the group of J3 values of 
us. Only one vertex is permitted to be chosen. 

(3) Nevertheless (9.2.1) would mean more for- 
brdden* points.) Constructvan’ axis Azk, perpendicular 
to si (Ooh eas) plane. The three axes AxK,, ae, 
and AZk, form a phase domain. A sphere in the 


Original K-domain (ref. Section 7.1) was described by 


2 2 aot Sy 
Ky + KS a oS ice i : (S39) 


Since our axes have been rescaled, (9.3.9) becomes 


(AxK,)“  (AyK,)*  (A2K,)* 
Ax Ay Az 
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which is the equation of an ellipsoid in the phase 
domaan” (Fag. 7). 
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PTSeOlmtnemsubbaces or the ellipsoid given by (9.3.10). 


and AyK are selected the value AzK 
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forbidden ellipses on the (AxK,, BES) plane.» “Only 
one point AS Mes bh xs Jar is allowed to be selected on 


such a family of ellipses. 
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Clarity thes, let us first consider the relation between 
CORO ee eos ao bands (9.302). 

Theorem 9.3 Provided AXK, 4/1 and AYE a are 
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rational numbers, and Ké - K 
PeLrrect souare, sthe tamily of ellipses, containing the 
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integers; p/Ax*, a/Ay*, d/Az* are rational numbers. 
Tiesonlyeesarvonal term 1S AZK4/1= ha fko-K 7K /M- 
Henceval etre ters tn (9.3.9) will not be summed up 
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For the second group of J values of K, namely 


3 3" selection will proceed in 


a .Samilam manner. Chesonly two differences are: 


the Ky for Ushi Spies 220); 
i] 
©3341 
is no longer arbitrary and should be equal to the Sy 
value selected previously in Section 9.3; and ee is 
Srbhebaryeput sheuldgnotsbesequalyto ry in order to 
avoid a singular R(J,)- 
ae ee = 
Ui Sheds" UUs eo) ee 
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values each, and they should be equal to Tyrloresen 


ie respectively. 


J1 

We shall proceed until we assign all Se mewn eve 
the rows. When this process is completed, we shall 
have J. different values of Si5 and Jy Cite rene 


values of Tha: Wernaverco far no reserictions for 
these values. But we understand that once we have 
chosen the Ja AxK, values and J5 AyK, values for r 


and s, all the i will be determined by our method 


just described. 


9.5 Stabilizing the Vandermonde Matrices 
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We know that the selection of AxK and AYR will 
eventually result ina series of matrices me (33) 4 

m = 1,2,6+2754551 which (see 8.4.3) are Vandermonde 
matrices belonging to Sequence C. To avoid a singular 
mz (33) being created, AxK,., and AYK 4 should in general 


not equalmasmul taplepoteie(Scegsectlone9. 3c). 
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We have so far not been concerned with the actual 


values of AxK 
; Sk 


AxK, and AyK a These actual yalues will be discussed 


and Ae hence the actual values of 


in Section,976-and,Chapter x1. 
We have not yet mentioned the choice of values ish 
and que although on some occaSions we said that Py and 


g,.. are both integers. In this section we shall describe 


u 
how to find the optimal Py and q, to guide our choice, 
To start the description, let us refer to Section 
4.2, in which it was indicated that Vandermonde matrices 
are known to be ill-conditioned and attention should be 
DesdalOmnstuahalit yOu salm: VS touse lect Py and qa 
such that me (33) can be stabilized to the best possible 
extent. 
Mhewanverston ol Sequence Cors given in (4.2.5). 
The matrix [D(n)D7*(n)] contains factors (x,-x,), 
Be (x 


“X3)y (xX5-x,) PoC bbe ae See OOECLOSed 


a ES 1 


LoRxoein value, then [D(n)D7*(n)] will have most of its 


Vs 
diagonal elements containing very small denominators, 
consequently the inversion of me (33) becomes unstable. 
instabidityawillealso, occur when fanyyone of the 

Xp rXore ee 1 X53 is equal to any of the others. To avoid 
this, we have to maximize the differences between any 
pair of them. 


In our case, the elements of me (33) are denoted 
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? nd >» We p 
vat : at 




















aia “sce sip) See nt ore : 

pel Iaeb. ‘y cela fice pee ee bos og) | 
Tk aoegnl) fine, 298 ba thie : — 

a dovtee) Loa ete, Any Spam omy sod ova om aes cy 
Ab git x ENOLSESSO Sie we Sadie “4? ts | 
og abt jal Otapante Med. vis wee 

42 rate sd (puibegs a& BagR-o2 wort - | 

300, = = \ lee oe ee wan7e oft 5 a) 

heme Stuy SRA bes ar cient 20 72 Fiekeyw az wS0d 
. ey oh Hos 3 pt oa ia one ste 

x ie 7s). «13 16 coda a2 tag ' 

icles ee i nes 4; En) “2 foye z 


_ 
: novip ait edges lie mpklebeel ae an 


rl oe a ule) Iga) Wil oes a gt ae 
sein sive: AES [ge] "aca ; tout ,allev nt o% of 
,cros gm mish Lheme eae ee eee ‘sdhemeta tanoysih: < 
eidtodans eumnagstl ight sh ptiVatesr xt ont 2¢) 
| bit. 29 “end pede. nee pace pris a". 
Dive ae Deis tath Rete 1. ¥ wire 


a a a tyes 





‘ah 
7 


118 


and call them ty and a0 Since 


th = exp(-i Az sobs ; (95Oe 2) 


ty Nema complexanumber. lO maximize the ditference 


between all the t's would mean to separate the t's 


WhEenmeduadleandlesewithan 2) (sceurigs. 9 and 10); 
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For simplicity, let us assume that we have 


selected Ky and we wish to see how we are going to 


select KorK3.-- PoOuMaiiime sLicestabimity.) in thas 
fated Daas OF Gee ue We can write ie SS} Toye Gkavcl Cha 
as q. We also assume that Ax = Ay = Az. (9.5.6) 
becomes 

Ie oslud « Me MeO A I 13 
Be el ae ae 8 ries ee oe 


Solving Lor cd, 


[a Oe pape? Bes 
“ ee (AyK.) 4m (p AxK, tmp lt/s 3 AZK 4/53) 
eye 
27 


(9793) 


We write p as Po and q as qo because the value obtained 
TowmanwovolLIMmalenialue,~moterthevactual sellectedivalue of p 
Saal Cops 

As a numerical example, we have the following 


experimental data: 
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AzK el Om racdtans < 
Fad 
Since AzK, >> AxK,, and AzK,, >> AyK 47 we can ignore 


the terms involving AxK 4 and ot Phe $0 oye, Leaving 


' 2 2 
~ + as we, 
qe t Vek /J4t - po - 1/955 (Onset) 
/ 2 
3 Joo = Pot ox 
Ye Lee 
an aa a = Waste | — as kee) 
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the above would be Pae 2 and = 2Peee nas. yatewds Fan 
optimal Ku amet riDuLcrony on the (AxK,, Does plane. 


Frome(Qeon means (rs 29.)c 
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AxK, 5 = AxK, 4+ AVGe SD, = ASN + 47 

(O05 FeES)) 
AxK, 3 = AxK, + 47 ; Gs = NS + 47 


AxK,  =AkK, + (n-1) 47, AyK Peg ae Bit : 


When choosing p and q in this way, we should not 
forget that the values of AxK, and Ayk, must represent 


a point on the ellipsoid. Thus, 
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TEPthis) requirementeis not met, p and q will be less 


than Py and do: 


OmoOmmEactOueMotrricese belonging loetie, Beoequence 


So far the values CyeeeeXsq7 Syeeee+Szo ,Corres- 


ponding to AxkK andy tor Ay 


x ((k,-1l)d WC) J+1) 


gre ie 
BOG ky = 1,2,..-5), ky = 1,2,++-J5, remain free. We 


can take advantage of these values to further expedite 


our computation. We know that r ie are theyr iest 


oat 
columnseor matrix R(J,) and Sye+++Sz5 are those of 
S(J.) ine ooo eee OU trast algoimiiom (3.4. 10), 
R(J,) has to be anverted and S(J5) has to be inverted 
and transposed. Perhaps a clever choice of these will 
make the inversions easy. 


One way to achieve this is to select (for another 


Wayjpsee Chapter x1) 


IC 1) dw SARS, (9.6.1) 
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This selection has avoided AxK, = -p7, thus avoided 

the AxK,-axis evimmetny 9 (rel. .9.5. 19)... Comparing "the 
apOve Witne (4.5.1), we can see that ry4 are identical 
tosthe first. column) of E(J,) ioetheereseduence: selias 


is followed immediately by the easy inversion 


= bier 
R (J,) = R(Jj) : (Oo) 


olielarlyewe shali= select 


A = 1(2k5-1)/d, : ComGre) 


yK eS 
y ((k, 1)J,53+1) 


This selection will avoid the ENA ean symmetry and 


Av) 
A) = si Gi (9.6.4) 


Therefore the solution to Pf = F in (8.4.10) becomes 


A A) 
— W 
Rect 55 Rect, f = R(J,) Rect>. E S(J5)/J5,J5- (9G. ) 
beteusespl te J. iitOl slacvo rs Jy and Jit 
Jz=J,Jc, Jas being positive integers . (ORO. 0) 
Then les 


P= (k3-1) Pp, ‘ k,=Int{ (u-1)mod (33) /5,}+1=1,2,..+755 * 


COR Oe) 


Gy = (Kg-1-75) 4, , k= (u-1)mod(J,)+1=1,2,...,54- (55 ot) 
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AXK = m(2k,-1)/d, tz 2(k3-l)pot 

C700) 
AYK 4, = T(2Ky-1)/J3, + 2(k,-1.75)a50 
ES es oh ie sY i (ky -l)JoJ4+ (ky-1)J,5,d, : Chae BoD) 


A numerical example will be given in the next 
chapter. Thetfollowing’is a summary of selection of 


parameters: 


Ji 1J5153 are determined by the physical model 
Lee 1 geass 

Jyrde are determined by (9.6.6). 

Ky rkorK3rkg Prd, are determined by (8.3.1), 


(One) bp leGry),) ja Or) ae where Po 
enascrna ree OO taInedm irom (O55) ang 
fe) 

u= 1,2,...5,5553. 
J41J9°33 Saonepemcal culaceOmurommos2 2) LOD 
vV=1,2,...5,J5553. 
AxK and AyK avemOotalnedme. roms (..6. o))s, 
cul yu 

AzkK pbacejih (US Si AW 0) )) 

ZU 


From the above we shall be able to determine r,s,t, 


[Sip oey li ign, Mees G 


The method of parameter selection helps to achieve 


the following: 
the 


(1) The matrix P created will comply with format 


Given by, (6.3). 4). 
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Singumeari ey Cine, oo OF me Matrices has been 
avoided. 

Seableliwcy ei eciTe ae Nace ces eS OPlLimi zed. 
MateLceseruanc —o be longetowchie Posequence, 
and are easy for inversion. 

The matrices me! hence matrix P,can be 
qutck ly sinvemted . 

We shall see in Chapter XII, that the best 
Cong Lemons nuMNbereEromeP is sobtainedstinough 
selecting phase angles in the ways described 


above. 
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Kee As COMPUTED EXAMPLE 


10.1 A Simplified Physical System 

Let J5= J5= J3= AS elicmscattering Object 15 
represented by 4x4x4 = 64 grid points, to each of which 
wemattacned a scatcered field source density foe Por the 
sake of easy comprehension, only 8 grid points are 
assigned a non-zero source field density £ = 2+i where 
— —Imand atethemolicm moines thersource densities are 
assumed nil. Fig. 11 shows the actual assumed positions 
of non-zero source fields as follows: 

2. BEOnm yee 27557 40,4 1, 507,51, 027,03. 
fo = CUO Rr eery 
O for—oiher values of v . 
We know we have continuous far-field data avail- 
abl eeeLGGmuSESeCOmloWwsWwesare Going = to choose a KL which 
leads us to select the data Fa at proper points and 


determine the matrix P for a quick solution of Pf =F. 


10.2 Selection of Phase Angles 


Let J3= Je =902..) From (9.6, LO). we nave 


u= k3 + 2(k,-1) a 4(k,-1) at (k,-1) (2x2x4) 


k < 2k, ate 4k, Ee 16k, Bie COE, Zierds) 


The k's will be computed from u according to (8.3.1), 
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1 ee for V=54, 35,46,47,,50°/51, 62,63. 
Lie =O for other values of v. 
FIG. || THE SCATTERER !S REPRESENTED BY EIGHT 


NONZERO SOURCE DENSITY SAMPLES 
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By assigning the k values as given in the Table 
above we have 64 values of Rue the projections of 
Wirehmontco: hie (AxK, ,AyK.) plane £orm a selection 
pattern as shown in Fig. 2 eelhewactldlavalles sO. Ky 
are calculated by computer and are listed in Appendix 


IMR Azk, is )e¢alewlated according to 


AzK = woe Viken — K - K 3 (LOR 2 R3)) 
ZU Oo a yu 


Computation es eae Spo can be done using (9.3.4) 
Bnan toes .one, “Computation of ty Wil pe pewLilmaccoroance 
With oalee) 9 eihese values form the matrices R4)e 
S(4) and mL CA). We store these matrices instead of 
SEO RING mE LOMSaVeCmsopace. 


In this example, we have 


r._ = expl-idxk,] = exp [-in (-4.25+4k,+0.5k,)] (LOR 2 64) 
= exp[-in(0.5k,-0.25)], ki= peo syn 

S,9 = exp[-iAykK, 1 = exp [-in (-2.25+2k,+0.5k,)] CO 20) 
= exp[-imn (0.5k,-0.25)], ky = Lee ia 

Therefore ee msi R(4) = $(4). The generated data 
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BiGui2 


PEO ERE CRIONSEAMRE RNG BORK Koy. Kaa 
AS PROJECTED ONTO THE (AxK,, AyKy) PLANE 


(The numbers attached to each point are the values of u) 


int a i a) 

4 6 a bs 
ride i ; Oy Pe ; 
| Be > ; 
o—-F a 

31 3! ri bh 






Ste | " mf Se} 
‘ae seid 
ey (et Ss i 
; 
a ee eee woe 
“ee f— a LES |. APSE, 
a eo 
met =v i 1. 





Tif kue Lied ,cwedl ‘tee 
(Cy iodeetiny sit ere ing” Ngee or 


y j 

oad - «= 

Tre F r. M @: 
. i : 


‘i ah 
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R(4) = =H (Ame LO 26) 
V2(-1+i)/2 -i Y¥2(1+i)/2 -1 


V/2(1+i)/2 a V2 (-1+i) /2 —1 


We can be sure that the inverse of R(4) or S(4) is 
R(4)/4 and the ordinary inversion algorithms are not 
needed. 3 

The og A) Maser cas erOr Na 2 ee, LO are@malso 
listed in Appendix III. Every me (4) sey. el (Ot Sieh ena es 
Its second column equals the square of its first 
Column, sand esethiisd "column equals the eine of its 
fixst column, etc: Hence “our fast inversion algorithm 
fomesequencenC awe apply. “Furthermore, *¥t “can be 
verified, from the numerical data, that ty, tor t3 and 
ty in every me (4) subtend big angles between them, so 
that anh Me) becomes well-conditioned as a result of 
adopting the optimized values for Py and que 

Pawisl ebemsamu lated, USING stheeR,. oo, as and the 
& Veieo sy (ONO) he) Se” Weave Wemlee kehslsyel 1S i eRILIN Gels) cyelra” Tenshers! 
for the system Pie= VF; 

Then we solve the system Pf = F by our fast 
algorithm to obtain £, and compare the obtained £ with 
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LOeo COMpULeT Programs: 


Our sconputang program 1S written in three parts: 

(A) A program to generate K-X components, R(4), 
Sta; pel 4s) and eo and) toe store the data 
in file DATA. 

(B) A program to solve Pf = F by fast algorithm. 
[ireecrogran 1s. Compl ledasand the ob ect 
program 1s) stored in file COMPUT. 

(C)mAn execution program to run COMPUT using data 
in DATA. 

Tiespregrameim Lub 1s attached ineAppendix EIT. 
Computer outputsparesalso,enclosed | §ilteistinteresting 
Comnotemthatacomsolve Pt = F using "a set OL given data 
requires only program (C), which needs 2.2 seconds 
CPU Ga Nem Os Jy= J5= J3= 4, The maximum error shown 
in the computed fT does not exceed 0.0006%, which should 


be acceptable to both scientific research and the optical 


industry. 
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MIS aVAReATLONS OF “THE, METHOD OF SOLUTION 


11.1 Factor Matrices Belonging to the Fourier Transform 


Sequence 


The algorithm described in previous chapters has 


been based on formula (8.4.10), namely, 
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AN ALTERNATIVE SELECTION PATTERN FOR K,...Kgq, 
USING THE FOURIER TRANSFORM SEQUENCE. 
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PART THREE 


ERROR ANALYSIS 






4 


rey AA 


Cleves 


AII. SEQUENCE CONDITION NUMBERS AND ERROR ANALYSIS 
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2. Whether we sforce Rand's into Sequence E, 
Cle McOmSeqlencemn,. the result willbe the same: 
IAS Se? ally NMSo si ale 

3. Sequences L and H have unity condition 
numbers. They should be ideal for inversion. 

4. The matrix T can easily become unstable. 
However, it can be optimized towards a Fourier 
Beant onocnieina tis! “aT ie whict CcaSsewica willl en Oy ay lini by 
condition number. 

Dei emaCCuLo CVO Ollig tas beclgonLEnmecai De 


expressed by an error bound 


latin  (llgaFll, | Udall, 


< (Cieree) 
[ltl 1. at es PRE | 


The er rac dound is n times the maximum relative error 
in the me Matrices plus n times the maximum relative 
errors in the data m° It can be interpreted as the 
N-condition number of P being almost a unity. This is 


the best P we can possibly create for the system of 


equations Pf = F, 


ret . 
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XLT = DLSsCuss LON 


This chapter summarizes our method of solution to 
inverse scattering problem and the mathematics developed 
in this thesis. The advantages and disadvantages of 
our approach as compared with other existing methods 


will be reviewed. 


TouleeGensless™ Reconstruction 


The conventional way of forming an image of an 
object through a system of lenses is known to be subject 
to the following conditions: 

(1) The object should be larger than the wavelength of 

EhemullLuniieeLigeradtat don. 

(2) The aberrations of the lenses are negligible. 
(3) A three-dimensional object results only in a two- 
dimensional image. 

The last point is considered as a major handicap 
of conventional imaging by lenses. 

Computational optics studies a wave field by 
numerical evaluation of field data. it needs no lens 
and forms no image, and is therefore not subject to the 
above conditions. With the help of mathematical analysis, 
the optical data obtained can be processed, and any 
particular information that enables us to characterize 


a wave field, for example, a spectrum of spatial 
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frequencies, can be extracted. The phase information, 

which no physical instruments of our time can detect 

directly, has to be sensed by an interferrometric 
method, like holography or a three-reference-beam 
method (see Schmidt-Weinmar 1973 and 1975). It is on 
this basis that computational optics leads to a new 
dimension in image reconstruction: the lensless recons- 
wre Siehe s(eValc 

Lensless reconstruction has some advantage over 

a conventional image forming system in that 

C1) the errors due to imperfect lenses are eliminated, 

(2) optical messages can be coded and compressed for 
transmission and can be decoded and restored at 
a remote distance, 

(3) noise can be filtered and distortions can be 
suppressed so that any optical information 
contained is enhanced, 

(4) a 3-D object can be reconstructed, 

(5) information can be extracted concerning waves of 
a particular type. 

This thesis gives an example for points (1), (4) 
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13.2 Reconstruction of Objects near the Resolution Limit 


Two successful approaches have been made previously 
to 3-D reconstruction. These were computer holography 
and, X=ray crystallography; both are lensless methods. 

Holographic techniques are limited to some appli- 
CAlionse. s lhiey  requine in particular that 
(ae) The spectral density Of the scatterer is bend-— 

limited to wienin a bandwidthwor Wes 27/8) os where 

ps is the vacuum wavelength of the illuminating 

light (see, for example, Wolf 1969). 

2) TiemsLecording: sadone: 10 theseresnet region of 
the scatterer, and over a plane, not a spherical, 
Sumrocer. 

)) Many sets of data are combined that are to be 
obtained at one fixed wavelength, and with 
different directions, of the light that illumi- 
nates the object. 

A holographic approach will thus not be applicable, 
if the object to be reconstructed has a microscopic or 
submicroscopic size. It is usually very difficult to 
take accurate recordings in the Fresnel region with a 
photographic film of ordinary thickness. 

X-ray crystallography TSeUuSse nul Only ebOrestLuCcture 
elucidation of crystals containing molecules of molecular 
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weight less than about 10°; Ehvced SeChewotetyOnOLlacemOL 


magnitude below that of many biologically crucial 
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molecules such as large nucleic-acid polymers (see 
Chapline and Wood: X-ray Lasers, Physics Today, June 
ih / hop) ie 

the approach developed in’ this thesis will be 
combined with the Enreen re rerenceeneon method to 
Compute, (romerar-rield= data, thiesspatteal distribution 
Of tnewscattercd field source density thet represents 
Spocabeerer., | Nemscatterca fieldrsource density as 
computed vine discrete Lorm. A density sample at some 
point in space represents a mean value over the source 
Censttyewi Citnmasolocklet inside thesscatterer.) This 
ion OtmeOomascumemiuliatstneeSscattererm=eitselt as a dic— 
crete object, but to use an approximate representation 
of the integral over a slowly varying continuous 
function by a finite number of: terms; each of these 
terms is a weighted average of the integrand in the 
interval belonging to the: particular term. From these 
samples of source density obtained, we may determine 
thewsSpacra leds tribution or Seniewareibn potential of the 
Objectpeive. ik (n* - te) LuLS OGbOVvidessa so lULdoOnmto 
the inverse scattering problem in an area where holo- 


graphic and X-ray methods cannot be applied. 


l3so bast inversion OL large Matrices 


A special feature of optical computation is the 


Manipulation of very large matrices, which may represent 


iff i =e: Poi i 
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either the image of an object or some distribution 

belonging to a propagated optical field, or some 

coded message from which an image or the distribution 
in space of an optical-wave field may be recovered. 

With the ever increasing demand for better image 

quality, ever larger matrices have been put in use. 

Conventional numerical methods are unable to handle 

such large matrices because of the exceptionally long 

CPU time and the large core storage required. The 

‘propagation of computing errors also presents a problem. 

Special techniques therefore need to be developed for 

large matrices. 

This thesis has demonstracted a new inversion 
technique outside the areas of banded, sparse, ortho- 
gonal, and ordinary Vandermonde matrices, to which so 
Many contributions have been made, The following 
results may be noted: 

(1) A number of fast methods for direct inversion of 
certain large matrices have been obtained. The 
rengemorematrices, introducedwinsGhaptexrs ia to 
IV, can be numerically mers eG as accurately as, 
and almost as fast as, orthogonal matrices. 

(2)) The inversion of certain composite matrices by 
parts has been shown. For example, a row-wise 


Kronecker matrix may be inverted, by multiplying 


ule 





em 
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(3) 


Ehemcolunns, of gam inverted factor matrix, with 
BiemoOtiemecenstituent matrices@minverted, 

A technique for solving a system of equations 
has been given based on rectangularization. An 
m-plex system (Section 6.3) of equations can be 
solved qurckiy, wif the data wand the unknown 
vectors are both rectangularized. 


The many examples given in Part I indicate that 


this is an interesting area, worthy of more intensive 


explorations gn sGueurel, 


13.4 


THiemMethodsot  s—-DeReCconstruction 


Our method of solution has been based on some 


physical assumptions: 


ey) 


(2) 


Ve) 


The object 1s weakly scattered and its refrac— 
ElLvVemindexsie varies slowly wLtninecne scatterer, 


14 ienweyer 


ieee nanny no 10 
inewsCaveecingobjJech 1s Ulluminaledeby a plane 
wave field that comes from a monochromatic and 
COHGLeNt Sounce. 

The data Ele are taken from a region far 
Way ee COMmUNeEOD|eCE Ts tice a] cum = ieee where 
XV iswarspacemVeccor ilethe far—imicideregmlon, sqlld 


Kis a space Vector fon any Sources point, that 
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(4) Data OE ee are available at any point 


1 

Kee (eae ee) ew che eee KE Ke) win 

eV. eZ, Z O x Vv 

the K-domain. 

Ours -Desecouseructiones jsimplified tovangunverse 
mapping of F, the far-field data, onto f, the scattered 
Pieldgsourcesdensity,a0r uns othe rewords,s tossolving 
asystemionvequations Prt = Fatfor: fee While ft remains 
inechiem(34,V,2) spaces dolain,nfsise redefined’ in the 
EEE NOS TEAS) domain. In view of K, being dependent 
on Ky. and Ke F may be considered as defined in the 
domain's projected plane a for convenience, 
The indices in the latter domain Ky Ko rkaree- are mapped 
ONnCOFy dmsui near index @Uj= 1:2)... nee Acmchersame. tame, 
the spatial indices J,1J9133 are also mapped onto a 
linear index v, v=1,2,...,n. When these indices are 
attached to F and f respectively, we have two Hilbert 
spaces withen-dimensional vectors F and) £, and the pro- 
pagatormel sMapSeimoncon teen lhe Kerncieot P) turns cute bo 
be separable and symmetric in the x-,y-,Z-directions, 
thus allowing the mapping of P-to be performed first in 
the x-direction, then the y- and finally the z-direction 
by thesoperators R, S and T repsectively. To reverse 
the mapping by P, it is then only necessary to reverse 
the three directional mappings, one at each time. 


The equi-spaced sampling in the (x,y,z) domain, 


as well as in the projected Wg Ea | plane, yields 
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Fredundanterectorsein matrex PS so that) P¥can”’ be 
decomposed to allow the application of fast matrix 
inversions. The selection pattern in the ey 
plane has J rectangular blocks where J, is the number 
of samples in space ire the z-direction. The distances 
between these blocks are equal to multiples of 27 in 
the directions of either AXK, ox avi the blocks are 
Ggiiedenticals each olock contains a mesh of J.J 


eZ 


points, where Jy and J5 are the number of sample points 
ines bacCeminethesx—-eand y-direct ions, 

Our fast algorithms begin by rearranging the data 
F, then post-multiplying them with Tee (see equations 
Span 7 ands 624.8) st his,is the inverse mapping inyz- 
direction, namely, we pick one point in a rectangular 
block, galdmtnes cocvresponding, points pinsevery other block, 
and operate the data F corresponding to these points 
with thesanverse—-map, operator Te As a result of this 
operation, we obtain a set of new data F", which are the 
F's already inverse-mapped in Zod recone Further 
inverse operations in the other two spatial directions 
can be done according to the formula Rect,,Rect,3f me 


i 


(Ro +) Rect Te Va lie, anverse operators vin stnelx—. and 


J 
y-directions are oo and re respectively. she rectangu-= 
larization of data before each operation, ensures that 


thespoints are aligned in the xight order. 
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Haye) 


After these three inverse-operations we shall 
obtain the required scattered-field source density f, 


and the equation Pf=F is solved. 


MeesOvo, alee culiios theimethod resented 


Our method of solution has been designed to yield 


an optimum result with regard to the following condi- 


UOIS 

(4) The solution of the linear system Pf=F is speeded 
up. 

G2) Realldsiomea Gem MaaEncesmrOtmwhichea method of fast 


PnVeEcsoi ONawomks. 


es) The Vandermonde matrices Ta? m= 1,2,..+,5,55, are 
Sheeler h. 

(4) Computation errors are minimized, 

(5)) Computer storage requirements are optimized. 


Our fast algorithm needs only J 455 (Iy+5,40 444.5 SG) 
complex arithmetic operations and approximately 
Ve gacth le complex words core ates This improve- 
ment allows us to determine, for example, a 80x80x20 
source density representation in 8 minutes on an IBM 
360/67 machine without virtual memory. Without such 


an algorithm, solution of a system o£ 128,000 equations 


would not be computationally feasible. 
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i? OE LoposedsDiveceL ous tor Future Research 


In view of the above discussions the following 
may be suggested for future research as a continuation 
ORE CSc nes hc. 

(1) To discover other fast-invertible matrices for 
the’ categories mentioned an Part I of the thesis 

(25) To develop other partial inversion techniques 
for composite matrices, 

(3 An-.experimental realization of the inverse 
scattering process to reconstruct the physical 
structure of a real object, 

(4) To,ycompute the optical field near and inside a 
scatterer given the magnitude and phase of the 
faratield, 

(5) Extension of the work towards reconstruction of 


aeSubmicroscopic Ob FEct . 


Vowel 
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APPENDIX I 


KRONECKER PRODUCT OF MATRICES 


DefinwetoneeA wKhronecker multiplication sis repre- 
sented by the operator @. Let A= ta, .}and B= Shes 
be two matrices of size (n,m) and (p,q) respectively. 


then theamulteipiication AS B yields 


iN GY BS fa, 5B} 


as a Kronecker product or Kronecker matrix of size 
(np,mg), which is expressible as a partitioned matrix 
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Decomposition of a Kronecker Matrix A Kronecker 
Matrix can be decomposed into an ordinary product of 
matrices which are sparse and of the same size as the 


original Kronecker matrix. For example, 
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These sparse matrices are usually referred to as Good 


matrices, and will have regularly repeated elements. 
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APPENDIX Lt 


DXKX(21 


DXKX (28) 


DX KX (29) = 
DX KX (30) = 
DXKX (31) = 
DX KX (32) = 
DX KX (33) = 


DX KX (34) 


.DXKX (35) = 
DXKX (36) = 
DX KX (37) = 
DX KX (38) = 


DXKX (39) 
DX KX (40) 


DX KX (41) = 
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0. 78540 
SIS h aE pe | 
0.78540 
1323519 
Ze IO 
14.92256 
Ze00d9 
14.92256 
3092099 
16.49335 
SaZoos 
16. 49335 
ree wate) 
18.06415 
5.49779 
18.06415 
0. 78540 
secu col 
0. 78540 
ec See 
Pee Meio) (Pe) 
14.92256 
Phy Biel) 4) 
14.92256 
Eyy SES! 
16.49335 
SeZoc9 
16. 49335 
Se Sila S 
18.06415 
5.49779 
78.06415 
0.78540 
3 ots Dail, 
0.78540 
Uses NTT 
2e550 19 
14.92256 
ZO Onis 
14.92256 
So SPAe es! 


THE 64 VALUES OF KY, SELECTED BY COMPUTER 


u 


= (AxK, AYK AzK 
DYKY( 1)= 0.78540 
DYKY (@2)= §0.78540 
DY KY (93) = 97.0685 
DYRY (94) = 97.06856 
DYKY( 5)= 0.78540 
DYRY{(36)= 0.78540 
DYRY (G7) =297.00655 
DYKY{ 8)= 7.06858 
LYKY( 9)= 0.78540 
DYKY(10)= 0.78540 
DYKY(11)= 7.06858 
DYKY(12)= 7.06858 
DYKY(13)= 0.78540 
DYKY(14)= 0.78540 
DYKY(45)= 77.06658 
DYKY(76)= 7.06858 
DYKY(17) =" 2.35619 
DYKY(18)= 2.35619 
DYKY(19)= 8.63938 
DYKY(20)= 8.63938 
DYNA 21) = 223 ONS 
DYKY (22) =" 92235619 
DYKY(23)= 8.63938 
DYKY(24)= 8.63938 
DYKY(25)= 2.35619 
DVKY (26)=e 2.5019 
DYKY (27)= 8.63932 
DYKY (28) = 8.65938 
DYKY(29)= 2.35619 
DYKY(30)= 2.35619 
DYKY(31)= 8.63938 
DYKY(32)= 8.63938 
DYKY(33)= 3.92699 
DYKY(34)= 3.92699 


DYKY(35)= 10.21017 
DYKY (36)= 10.21017 


DYKY (37) 
DYKY (38) 
DYKY (39) 
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DYKY (40)= 10.21017 
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DZKZ (15) 
DZKZ (16 
DZKZ (17 
DZKZ (18 
DZKZ (19 
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DZKZ (21 
DEEZ (22 
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DZKZ (25 
DZKZ (26 
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DXKX (42) 
DX KX (43) 
DX KX (4U) 
DX KX (45) 
DX KX (46) 
DX KX (47) 
DX KX (48) 
DX KX (49) 


DX KX (52) 
DX KX (53) 
DX KX (54) 
DX KX (55) 
DX KX (56) 
DX KX (57) 
DX KX (58) 
DXKX (59) 
DX KX (60) 
DX KX (61) 
DX KX (62) 
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DX EX (G3)= 


DX KX (64) 


— 
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16.49335 
Shope eA eve hs) 
16.49335 
5.49779 
18.06415 
5e49779 
18.06415 
0.78540 
er) Wague 
0.78540 
Ae eS fal 
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14.92256 
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14.92256 
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16. 49335 
S620 99 
16. 49335 
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18, 06415 
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DYKY (42) = 
DYKY (43)= 
DYKY (44) = 
DYKY (45) = 
DYKY (46)= 
DYKY (47) = 
DYKY (48) = 
DYKY (49) = 
DYKY (50) = 
DYKY(51)= 
DYKY (52)= 
DYKY(53)= 
DYKY (54)= 
DYKY (55) = 
DYKY (56)= 
BYRN (57) = 
DYKY (58) = 
DYKY (59)= 
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DYKY (61) = 
DYKY (62)= 
DYKY (63) = 
DYKY (64) = 
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67.68008 
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67.64006 
ey Shei chewarelt | 
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223.1085 
59 .89409 
61.-43959 
58.98080 
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ISN 
Low 
ISN 
ISN 
ISN 
HSN 
ISN 
ISN 
ISN 
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ISN 
ISN 
ISN 
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A. 


0002 
0003 
0004 
0005 
0006 
0007 
0008 
0c09 


0010 
0011 
0012 
0013 
0074 
0015 
0016 
0017 
0018 
0019 
0020 


0021 
0022 
0023 
0024 


DEPEND TX) Ane. 


AY Computer) Program to Generate Matrices R, 65,,.T, 


Ano a ca Eeandeotores thes Outputs im) Fide DATA: 


Cc 
Cc 


20 


@QIQE@Qr@rG 


150 


170 
180 


50 


*OPTICNS IN EFFECT* 


*OPTICNS IN EFFECT* 


*STATISTICS* 


*STATISTLCS* 


NO 


4bsKOGK PROGRAM TO GENERATE DATA #**#%%%% 


COMPLEX R(4,4),7(16,4,4) 
COMPLEX CMPLX,F(16), BIGF(64) 
COMMON NX,NY,NZ,TWOPT 

READ (5,20) NX,NY,NZ,OK 
FORMAT (312,F5.2) 

NXY=NX*NY 

NXYZ=NXY*NZ 

TWOPI = 6.2831853 


COMPUTER SIMULATED DATA 

ACTUAL PROBLEM BIGF WILL BE MEASURED. 
MATRICES R,S,@ WILL BE COMPUTER GENERATED. 
BIGF IS NCW COMPUTER GENERATED. 

CALL GENR (R, NX) 

WRITE (6,150) 

FORMAT ("1") 

HOminLONI=1;6 

WRITE (6,180) 

FORMAT (t=') 

CALL GENT (T,NXY,OK) 

WRITE (6,150) 

CALL GENBF (BIGF,R,T) 

WRITE(7,50) ((R(J1,d) ,J=1,NX) ,J1=1, NX) 
WRITE (7,50) (((T(K1,K2,K3) ,K3=1,NZ) ,K2=1, 
NZjPea—1), XY) 

WRITE(7, 50) (BIGF(1) ,1=1,NXYZ) 

FORMAT (i (E1 Sent Ed Sele) 

STOP 

END 


NAME= MAIN,OPT=01,LINECNT=59, 


SOURCE, EBCDIC, NOLIST, DECK, NOLOAD ,NOMAP 


SOURCE STATEMENTS = 23 ,PROGRAN Sisk = 


DIAGNOSTICS GENERATED 


ea END OF COMPILATICN *¥4#4*** 


COMPILER STATISTICS: ELAPSED TIME 2.620 SEC. 
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MTS FORTRAN H (OS REL 21.7) (CT204) 


COMPILER OPTIONS - NAME= MAIN, OPT=01,LINECNT=59, 
SOURCE, EBCDIC, NOLIST, DECK, NOLOAD, NOMAP 


Cc 
ISN 0002 SUBROUTINE GENR(R,NN) 
Cc GENERATION OF MATRIX R & S, R=S IN THIS CASE. 
ISN 0003 COMPLEX R(NN,NN),CMPLX 
ISN 0004 A=3.141593 /NN 
ISN 0005 DO-10.J31=4,NN 
ISN 0006 JJ= (2*J1-1) 
ISN 0007 Doms Omd— 1 NN 
ISN 0008 COSA = COS(¥Jd*J) 
TSN 0009 IF (CCOSA.LT.-0. 999999) COSA=-1, 
Isn 0011 IF (COSA.GT.0.999999) COSA=1, 
ISN 0013 LF ((COSA.LT.0.000002) AND. (COSA.ST.-0.000002 
N\NMCOSH=0-.0000 
ISN 0015 SINA =-~SIN(A¥JJd*J) 
ISN 0016 IF (SINA.GT. 0.999999) SINA=1. 
ISN 0018 LF ((SINAe LT. 0.000002) «ANDe (SINA.GT.-0.000002 
-)) SINA=0.0000 
ISN 0020 IF (SINA.LT.-0.$99999) SINA=-1- 
ISN 0022 10 R(J1,d). = CMPLX (COSA,SINA) 
ISN 0023 WRITE (6,150) 
TSN 0024 150 FORMAT ('1') 
ISN 0025 Domig0 £-1,4 
ISN 0026 170 WRITE (6,180) 
ISN 0027 180 FORMAT ('-') 
ISN 0028 WRITE (6,30) ((R(J1,5) ,J=1,NN) ,-d1=7, NN) 
ISN 0029 30 FORMAT ("'OR=",4(4(F7.4,F704,1X) /e' ')) 
ISN 0030 RETURN 
ISN 0031 END 
*OPTIONS IN EFFECT* NAME= MAIN, OPT=01,LINECNT=59, 
*XOPTICNS IN EFFECT* SOURCE, EBCDIC, NOLIST, DECK, NOLOAD , NOMAP 
* STATISTICS* SOURCE STATEMENTS = 30 ,PROGRAM SIZE = 1090 


*KSTATISTICS* NO DIAGNOSTICS GENERATED 
4444%xX 0 ENDSCY COMPILATIONS *4*F%~ 


COMPILER STATISTICS: ELAPSED TIME 16112 GOBOEC. 
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MTS FORTRAN H (OS REL 21.7) (CT204) 


COMEILER OPTICNS - NAHE= MAIN, OPT=01,LINECNT=59, 
SOURCE, EBCDIC, NOLIST, DECK, NCLOAD, NOMAP 


ISN 0002 SUBROUTINE GENT(T,NXY,OK) 
C GENERATION OF MATRIX 7. 
ISN 0003 COMMON NX,NY,NZ,TWOPI 
ISN 0004 COMPLEX T(NXY,NZ,NZ), CMPLX 
ISN 0005 OKSQ=OK*OK 
ISN 0006 NZD OM = ENA 2 
ISN 0007 A=3.141593 /NX*2. 
ISN 0008 B=2.141593 /NX*2. 
Cc FOR SAKF CF SIMPLICITY DXKX IS WRITTEN AS XK, 
E WEST ING SEX- PAGS NS) “ASS 
ISN C009 YK1=-3.141593 /NY-6. 2831853 
ISN 0010 XKO=-3. 141593 /NX-12.56637 
ISN 0011 N=0 
ISN 0012 NNT = 0 
ISN 0013 D CHI Omk2 1 NX 
ISN 0014 YK1=YK14B 
TSN 0015 XK1=XKO 
ISN 0016 DOMTA0K1= 1, NX 
ISN 0017 XK1=XKK1+A 
ISN 0018 NNT= NNT +1 
ISN 0019 NT=0 
ISN 0020 YK=YK1 
ISN 0021 DOM 1OmKo= 12 
ISN 0022 YK=YK+6.283185 
ISN 0023 XK=XK1 
ISN 0024 DO 110 K3=1,NZD2 
ISN 0025 XK=XK+12. 56637 
ISN 0026 NT=NT+1 
ISN 0027 N=N+1 
Cc ASSUMING DX=DY=DZ=1 MICRON. 
Cc ACTUAL FORMULA IS ***  DZKZ = 
C DZ* (SORT (OKSQ= (DXKX/DX) **2= (DYKY/DY) **2) ) 
ISN 0028 ZK=SORT (OKSQ=XK*¥XK-YK*YK) 
ISN 0029 WHITTEN G20 )meNEMe, CK NMA Ne, os 
ISN 0030 120 FORMAT({* DLR oso LR GG 
tips pina vaca Cates Uy aU ahs Se 
ISN 0031 DO 130 K=1,NZ 
ISN 0032 COSA = COS(ZK*K) 
ISN 0033 SINA = -SIN(ZK*K) 
ISN 6034 430 T(NNT, NT,K) = CMPLX ( COSA, SINA ) 
ISN 0035 410 CONTINUE 
ISN 0036 — WRITE (6,150) 
ISN 0037 450 FORMAT ('1') 
ISN 0038 DCm Omiaiy 
ISN 0039 170 WRITE (6,180) 
Isn 0040 180 FORMAT ('-') 
ISN 0041 HOM OONK 1-116 


ISN 0042 460 WRITE(6,140) ((T(K1,K2,K3) ,K3=1,NZ) ,K2=1,NZ) 


Ok 
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ISN 
ISN 
ISN 
ISN 
ISN 


0043 
0044 
0045 
CO46 
0047 


440 


200 


FORMAT (*OT=',4 (4 (F7e4,F7e4, 1X) /,' 
WRITE (6,150) 

HomzZo0mI= 11,2 

WRITE (6,180) 

Dom190) K1= 7, NXY 


UP) 
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ISN 0048 190 WRITE(6,140) ((T(K1,K2,K3) ,K3=1,NZ) ,K2=1,NZ) 
ISN 0049 | RETURN 
ISN 0050 END 
*XOPTIONS IN EEFECT* NAME= MAIN,OPT=01,LINECNT=59, 
*OPTICNS IN EFFECT* SOURCE, EBCDIC, NOLIST, DECK, NOLOAD, NOMAP 
* STATISTICS* SOURCE STATEMENTS = 49 ,PROGRAM SIZE = 1792 


RSE ALtotrCos | NOY DLAGNOSTACS GENERAZDED 
AXE END GF COMPILATION 444%%% 


COMBI TRH  SpALL ST LCs: ELAPSED TIME 4,626 SEC. 
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MTS FORTRAN 4H 


ISN 


ISN 
ISN 
ISN 
ISN 
ISN 
ISN 
ISN 
ISN 
ISN 
ISN 
ISN 
ISN 
ISN 
ISN 
ISN 
bes N 
ISN 
ISN 
ISN 
ISN 
ISN 


ISN 
ISN 
ISN 
ISN 
IsN 
ISN 
ISN 


0002 


0003 
0004 
0005 
0006 
0007 
0008 
0609 
0010 
0011 
0012 
0013 
0014 
0015 
0016 
0017 
0018 
0019 
0020 
0021 
0022 
0023 


0024 
0026 
0027 
0028 
0029 
0030 
0031 


a 


(OSmRELM et o1) (C'T204) 


COMPILER OPTIONS - NANE= MAIN, CPT=01,LINECNT=59, 


OIGr@i@r@uG 


— 


430 


Lu0 
450 
451 


*OPTIONS IN EFFECT* 


*OETICNS IN EFFECT* 


*STATISTICS* 


SOURCE, EBCDIC, NOLIST,DECK, NCLOAD, NOMAP 


SUBROUTINE GENBF(EIGF,R,1) 

A ROUTINE TO SIMULATE THE DATA BIGF 

FOR SAKE OF SIMPLICITY, A SCATTERING OBJECT 
IS REPRESENTED EY SCATTERING POTENTIALS AT 
POINTSe@iz, 19S) e2p 47s, wert e 2.874), 
(S70 eas hoe (ro7 2) reset) 6 PHERE 
SMALLF = 2+, I**2=-1, OTHERWISE SMALLF = 0. 
COMMON NX,NY,NZ,TWOPI 

DIMENSION L1(8) ,L2(8) ,L3 (8) 

COMPLEX BIGF (64) ,SMALLF 

COMELEX R(NX,NX) ,CMPLX,17(16,NZ,NZ) 

Dome Ie 

WRITE (6,2) 

FORMAT ('-') 

NXYZ = NX¥NY*NZ 

N=0 

SHRLEPE-SOMDLN (col) 

DOmd00mi-4, 6 

rao ERG, Tihs Aieesy saved, 

FORMAT (311) 

N=N+1 

BIGF(N) = CMPLX (0.,0.) 

NNT =(N-1) /NZ+1 


NR = MOD((N-1)/NZ,NX) + 1 
NS = (N-1)/(NX*¥NZ) + 1 
NT = MOD(N-1,NZ) +1 


DO 430 1=1,8 

BIGF (N) =BIGF (N) +R (NR,L1 (Z)) *R (NS, L2 (I)) ¥T( 
NNT,NT,1L3(I)) *SNALLF 

IF (N.GE.NXYZ) GO TO 440 

GO TO 420 2 

DO 450 I=1,NXYZ 

WRITE(6,451) I,BIGF (TI) 

FORMAT( 10X,% BIGF(',12,°%)=',F15.7,F15<7/) 
RETURN 

END 


NAME= MAIN, OPT=01,LINECNT=59, 


SOURCE, EBCDIC, NOLIST, DECK, NOLOAD, NOMAP 


SOURCE STATEMENTS = | 30 ,PROGRAM SIZE = 1242 


*STATISTICS* NO DIAGNOSTICS GENERATED 


444% END OF COMPILATION ****%% 


COMPILER 


STATISTICS: ELAPSED TIME list 3G eS bG.s 
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The Matrix R(4) or S(4) Generated by Computer 


20 -1.0000 -0.7071-0.7071 ~1.0000 0.0 
20 1-0000 0.7071-0.7071 -1.0000 0.0 
0 =1.0000 0.70071 0.7074 =41.0000 0.0 
«9 4.0000 -0.7071 0.7071 -1.0000 0.0 


R= 0.7071-0.7071 0 
— (Vee OF. 1-0. 07 10 
=O. 08 0 00) 9.0 

Oe Oui 0.7074 =O 
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T= 0.9967 
0.0579 
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-0. 4u864 


T= 0.9928 
-0.3006 
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The 16 Matrices T(4) Generated by Computer 
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The 64 Values of Be Generated 


BiGF{ 1)= -6.2256298 
BIGF( 2)= -0.5613050 
BIC (ns) = -2.2156420 
EFIGF( 4)= - 3. 7373371 
BIGE (5) le NOSE 
BIGF( 6)= 0.9054432 
EPGE(e7)= Qai614275 
Bick (ao)= Oe a 227 
FIGF( 9)= 0.6133854 
PBIGF (10) = 1.0779219 
BIGF(11)= 2.5166521 
Bac (a2) — 0.2893295_ 
BILGE (f3)i— 0.8399563 
BIGF (14)= -1.2671318 
PG eo) — =5.6654854 
BIGF(16)= 0. 0395203 
FIGE (A7)= ~10. 3667536 
BIGF (18)= 5.5610733 
BIGF (19) = -6.4317226 
BIGF (20)= -7.6865349 


BIGF(21)= 0.6263094 


by Computer 


1.1064100 
4.6524277 
-5.7637615 
0.1444349 
~-1.9204302 
-71.3089981 
123455830 
0.5006943 
-2.5383072 
-0.3698998 
Oei22isda 
0.5678787 
6.2070236 
-0.6242361 
(1.8798981 
-0.1162157 
11.1930933 


9.6120625 


13.0047417 
0.8837890 
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BIGF (22)= -0.1444349 
BIGF(23)= Sq eevee 
BIGF (24)= 2.0007019 
EIGF(25)= -2.5111609 
BIGF (26)= 1.4820700 
BIGF (27) = Smo 62 542 
BIGF (28)= 0.4620686 
EIGF (29) = 10.6225891 
PIGF (30)= -3.0517349 
BIGF(31)= -12.6921244 
BIGF (32)= -0. 33763414 
BIGF (33)= -3.5750589 
BIGF (34)= 9.1402493 
EIGF(35)= -5.5165701 
EIGF (36)= -6.0702972 
BIG? (37) = ~2.5111656 
BIGF (38) = -1.81€67 04 
BIGF (39) = 4.9420080 
BIGF (40)= 1.2671318 
BIGF (41)= -4.9829588 
EBIGF (42)= 0.2333988 
BIGF (43) = 5.5622482 
BIGF (44) = 0.0542316 
BIGF (45) = 14,3868771 
BIGF (46) = -1.9313088 
PIGF (47) = -12.1498728 


BIGF (48) = -1.0951471 


= Se oon 
2.6768169 
0.8067845 
=e OOS ted 
= 5 WISTS) 
=O. 2490253 
0.7744809 
1056993713 
0.5681009 
Died Oeil ah 
1e23648845 
14.7943602 
Deo 12/983 
| Ac EAT PASS) 
0.2750196 
S18 Benes) 
93 @ O22) 0 
267957745 
0.6242361 
She ees 
SE MOTE ST 
0.0474453 
0.1101189 
4.1374674 
125677176 
4,5382385 


229674492 


LS 








bbibs hast t i 


oe =a 1 svia ‘ag 


BIGF(49)= 
BIGF (50)= 
BIGF(51)= 
BIGF (52) = 
BIGF (53) = 
_ BIGF (54) = 
EIGF(55)= 
EIGF (56)= 
BIGF(57)= 
BIGF (58) = 
BIGF (59) = 
BIGF (60) = 
BIGF(61) = 
BIGE (60) = 
EIGF (63)= 


BIGF (64) = 


028399525 
4.0722971 
-1.7449089 
-1.6398525 
-7.8225460 
me 039/530 
1.6471977 
0.7654558 
-2. 4683943 
= Ora 1 5935.0 5 
2.0643368 
-0.0742826 
6.09712542 
-0.3284135 
~4.7900763 


~1.2231808 


6.2070198 
0.7445297 
~5.1794376 
-0.2827377 
-1.8357210 
-0.8253107 
1.4098988 
0.13377900 
~ 0. 7098782 
-0.8149223 
024247271 
-0.3194571 
-0.5498953 


6.7790203 


1.8246908 
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AeCourcemerogramato solve, Pip= Fr 


Dig S epeog Gam srSato le Sena Rae under *FORTG and the. 


object program will be were igh ARAL Sy YOKOVMID) Jes Reyes ed eh eve eablisyers 


ISN 
ISN 
ton 
ISN 
ISN 
JAY 
ISN 
ISN 
ISN 
ISN 
ISN 
ISN 
ISN 


ISN 


ISN 
ISN 
ISN 
ISN 
ISN 
ISN 
Eon 
ISN 
ISN 
ISN 


ISN 
ISN 


ISN 
ISN 
IsN 
ISN 
ISN 
ISN 
ISN 
“LON 
ISN 
ISN 
TSN 


0002 
0003 
0004 
0005 
0006 
0007 
0908 
0009 
00106 
0011 
0012 
heh Bs) 
0014 


0015 


0016 
0017 
0018 
0019 
0020 
0021 
0022 
0023 
0024 
0025 


0026 
0027 


0028 
0029 
0030 
0031 
0032 
0033 
0034 
0035 
0036 
0037 
0038 


C 


50 


310 


350 


150 


eek KKKKEK MATN PROGRAM 2K KK KE SK Ae 36 He 
COMPLEX R(4,4),7T(16,4,4) -U(16,4) -V (4,4) 
COMPLEX CMPLX,F (16), BIGF (64) 

COMPLEX M (4,4) ,Q(4) , SMALLF (64) 

COMMON NX,NY,NZ,THWOPI 

NX=4 

Ny=4 

NZ=4 

NXY=NX*NY 

NXYZ=NXY*NZ 

TWOPI = 6.2831853 

READ (7,50) ((R(d1,d) -J=1,NX) ,51=1, NX) 
FORMAT (4 (F15.7,F15.7) ) 

READ (7,50) (((T(K1,K2,K3) ,K3=1,NZ) -X2=1, 
NZ) ,K1=1, NXY) 

READ (7, 50) (BIGF(I) ,I=1,NXY2Z) 
SOLUTION OF PF=BIGF FOR F, USING FAST ALGORM. 
AS A FIRST STEP, MULTIPLY BIGF WITH TM1. 
NF =0 

DOmSs0mI—1 ,NXY 

CALL TRANSF(BIGF, NF, F,NZ) 

fe) hc) Beare 

Dome ONKad ANZ 

WALLS) = Weep TAS) 

CONTINUE | 

CALL INVV (V,NZ,M) 

CALL MATMLT (M,F,0 ,NZ,NZ,1) 

DOs Ond—1 NZ 

TRANSPOSING U 

U(i,dJ) = Q() 

CONTINUE 

SOLVING SYSTEM COLUMN BY COLUMN 

NF=0 

POMS 50m 1 NZ 

DO bs GON I= 1 Ney 

yen) Gees 

CALL LAMSCL(F,NXY,R, SMALLF, NF) 
CONTINUE 

WRITE (6,150) 

FORMAT ('1') 

CALL PROUT (SMALLF ,NXYZ,1) 

STOP 

END 





7ORPL ONS ON ME RP EGT* Nat E— MAIN, OPT=01,LINECNT=59, 
*OPTUONS IN EERECT= SOURCE, EBCDIC, NOLIST,DECK, NOLOAD, NOMAP 
SPA tow Cox SOURCE STATEMENTS = 5 /e, PROGRAM SIZE = 


eS eA Pollo NOM DERGNOsSTICS GENERATED 


oe eee ND Ore COMPS LAT EC Newark 


COMPTLER STATISTICS. ELAPSED Taue Boo 0ns ECs 


1 oh 


Dee 


A) ‘ai ri. poe 


ciroe te. 
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MIS FORTRAN H (OS REL 21.7) (CT204) 


COMPILER OPTIONS - NAME= MAIN, OPT=01,LINECNT=59, 
SOURCE, EBCDIC, NOLIST, DECK, NOLOAD, NOMAP 


Cc 
ISN 0002 SUBROUTINE TRANSF (BIGF,NF,F,NZ) 
Cc A ROUTINE TO TRANSFER THE ELEMENTS OF BIGF, 
Cc FROM (NF+1)TH TO (NF+NZ)TH, TO THE F(NZ). 
ISN 0003 COMPLEX F(16),BIGF(64) 
ISN 0004 ON 10) Lele NG, 
ISN 0005 NF = NFt1 
ISN 0006 10 F(I) = BIGF(NF) 
ISN 0007 RETURN 
ISN 0008 END 
*OPTIONS IN EFFECT* NAME= MAIN, OPT=01,LINECNT=59, 
*xOPTIONS IN EFFECT* SOURCE, EBCDIC, NOLIST, DECK, NOLOAD, NOMAP 
*xSTATISTICS* SOURCE STATEMENTS = 7, BROGRAME Siece = 356 


*STATISTICS* NO DIAGNOSTICS GENERATED 
axe END OF COMPILATION **¥*%* 


COMPILER STATISTICS: ELAPSED TIME UKE SRS 
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MTS FORTRAN H (OS REL 21-7) (CT204) 


COMPILER OPTIONS - NAME= MAIN, OPT=01,LINECNT=59, 


C 
ISN 0002 

C 

G 
ISN 0003 
ISN 0004 
ISN 0005 

ISN 0006 50 
ISN 0007 
ISN 0008 


*OPTICNS IN EFFECT 


*OPTIONS IN EFFECT* 


SOURCE, EBCDIC, NOLIST, DECK, NOLOAD, NOMAP 


SUBROUTINE INVE(E,N, Q) 

A ROUTINE TO INVERT A SPECIAL MATRIX E. 
HERE Q IS EQUAL TO E ROTATED CLOCKWISE. 
COMPLEX E(N,N) ,Q (N,N) 

DO)50)J3=1,.N 

DOP 50nd=1 9 

Q (3,03) =E(N -J3+1,d) /N 

RETURN 

END 


NAME= MAIN, OPT=01,LINECNT=59, 


SOURCE, EBCDIC, NOLIST, DECK, NOLOAD, NOMAP 


cate i Rm Roy Oho SOGURCE SSTADENEND Ss — H/ plMiselercaNG SSyaiAei = 562 


*STATISTICS* NO DIAGNOSTICS GENERATED 


HKKKKK END OF COMPILATION *%*¥%* 


COMPILER STATISTICS: ELArOL De LINE A E2st SIRE 
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MTS FORTRAN H (OS REL 21.7) (CT204) 


ISN 


ISN 
ISN 
ISN 
ISN 
ISN 
ISN 
ISN 
ISN 
ISN 


0002 


0003 
0004 
0005 
0006 
0007 
C008 
0009 
0010 
0071 


OP TON Seen 


SOP PEON SUNSET RECT 


-STAL to TlCs* 


COMPILER 


10 


EFFECT* 


Mh elo LCS ae eNO 


OPTIONS - NAME= MAIN, OPT=01, LINECNT=59, 


SOURCE, EBCDIC, NOLIST, DECK, NOLOAD, NOMAP 


SUBROUTINE MATMLT(A,B,C,M,N,L) 
THE SUBROUTINE WILL COMPUTE MATRIX C(M,L) AS 
A PRODUCT OF A(M,N) AND B(N,L). 
COMPLEX  A(M,N),B(N,1) ,C(M,L) 
COMPLEX CMPLX 

D019 l=, M 

DONO) d=1 

G (iy =CHPE (0s,0'e) 

DOM1ONK=15N 

C (I,J) =C (I,J) +A (I,K) ¥B(K,J) 
RETURN 

END 


NAME= MAIN,OPT=01,LINECNT=59, 


SOURCE, EBCDIC, NOLIST, DECK, NOLOAD, NOMAP 


SOURCE STATENENTS = 10 ,PROGRAM SIZE = 645 


DIAGNOSTICS GENERATED 


AkeKKK END CF COMPILATION **#4*** 


COMPILER STATE ST iCSs ELAFSED TIME le dch SHES 


had 










‘ereer a ae rievery 


pt, tee eee 
Ae i} herds Bo 
. ae) eat ‘tek 

SiS. (2 .<re 





os) 88 Ry oF 
a 
. xr " , 
Sead ie cr A et Kaw = | 
ra ey | “ a i. * 2 « ’ 7 
1avod, LOOM t98d CAL JON See LSRURE 7 rane 
: ee “% =36 7 * =| « + ~ 


esivivaspet:, or) alate heli 
Teh ca aor OupAT bu” 
cemeee en tRttbaon va 
—aT oo) (eeenis  <eorrexpare # 





NTS FORTRAN H (OS REL 21.7) (CT204) 


COMPILER OPTIONS ~ NAME= MAIN, OPT=01,LINECNT=59, 


18 


SOURCE, EBCDIC, NOLIST, DECK, NOLOAD, NOMAP 


Cc 

ISN 0002 SUBROUTINE PROUT(C,M,L) | 

| fe THE SUBROUTINE WILL OUTPUT THE MATRIX C 

ISN 0003 COMPLEX  C(M,L) 

ISN 0004 DOmieei= t,o 

ISN 0005 1 WRITE (6, 2) 

ISN 0006 2 FORMAT ('-') 

ISN 0007 DO Ope =i 

ISN 0008 10 PRIN THoOr ad 3(C(2,.)) 7021 ,L) 

ISN 0009 20 FORMAT (10X,° P(e =A 15 ee 15.7 /} 

ISN 0010 RETURN 

ISN 0071 END 
*xOPTICNS IN EFFECT* NAME= MAIN,OPT=01,LINECNT=59, 
*xOPTIONS IN EFFECT SOURCE, EBCDIC, NOLIST, DECK, NOLOAD ,NOMAP 
*xSTATISTICS* SOURCE STATEMENTS = 10 ,PROGRAM SIZE = 


woLAttot LCS* ~NO DIAGNOSTICS GENERATED 
HARE END OF COMPILATION *%*%%* 


COMPVEER SE SLAP StL Cs: ELAPSED TIME 5 AS SAGA 


434 


eet 










1. .. Meas Gaz 
+o GeE. ,45 S14, oh, eo epee 
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Mis FORTRAN TH (0S TREL (271.7) (CT204) 


COMPILER OPTIONS = NAME= MAIN, OPT=01,LINECNT=59, 


ISN 0002 


ISN 0003 
ISN 0004 
ISN 0005 
ISN 0006 
ISN 0007 
ISN 0008 
ao 30009 
ISN 0010 
ISN 0011 
ISN 0072 
HS Na 0203 
ISN 0014 10 
ISN 0015 

ISN 0016 


Or EONS IN EPPECT* 


A UbONSee Nebr ECDs 


SOURCE, EBCDIC, NOLIST, DECK, NCLOAD, NOMAP 


SUBROUTINE LAMSOL (Y, NXY,R,X,N) 

A ROUTINE TO COMPUTE F=(RM1) (U) (SN1T) 
COMMON NX,NY,NZ, TWOPI 

COMPLEX R(4,4),0(4,4),SM1T (4,4) 
COMPLEX ¥ (NXY),H(4,4) ,F (4,4) ,0 (4,4) ,X (64) 
CALL INVE (R,NX,Q) 

CALL ARANGE(Y,NXY,U, NX, NY) 

CALL INTRAN (R, NX, SM1T) 

CAP LEN Atta UUG otis, MONX NY, NY) 

CALL MATHLT(O,M,F,NX,NX,NY) 

DO 10 J=1,NY 

DOM Oer= 17 NX 

N=HN+1 

X (N) =F (1, J) 

RETURN 

END 


NAME= MAIN,OPT=01,LINECNT=59, 


SOURCE, EBCDIC, NOLIST, DECK, NOCLOAD ,NOMAP 


Me neo CO SOURCE STATEMENTS = 15 ,PROGRAH SIZE = dee 


*STATISTICS*® NO DIAGNOSTICS GENERATED 


Soot Ne Ore COMPACT AT EC Nex sa 


GONEP@L ERS SIA Leo tla. ELAPSED TIME -806 SEC. 


3: 
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MTS FORTRAN H (OS REL 21.7) (CT204) 


COMPILER OPTIONS - NAME= MAIN, OPT=01,LINECNT=59, 
SOURCE, EBCDIC, NOLIST, DECK, NCLOAD, NOMAP 


Cc 
ISN 0002 SUBROUTINE ARANGE (X,N,XPRIME,NBYJ,J) 
: Cc A ROUTINE TO REARRANGE A (N,1) MATRIX TO 
Cc FORM A(N/J,J) MATRIX. 

ISN 0003 COMPLEX X(N) ,XPRIME(NBYJ,J) 

ISN 0004 i=0% 

ISN 0005 DO 10 NN=1,d 

ISN 0006 DO 10 MM=1,NBYJ 

ISN 0007 I=1+1 

ISN 0008 10 XPRIME (MM,NN) =X (I) 

ISN 00069 RETURN 

ISN 0010 END 
*OPTICNS IN EFFECT NAME= MAIN, OPT=01,LINECNT=59, 
*OPTIONS IN EFFECT* SOURCE, EBCDIC, NOLIST, DECK, NOLOAD, NOMAP 
* STATISTICS* SOURCE STATEMENTS = 9 ,PROGRAM SIZE = 440 


*STATISTICS* NO DIAGNOSTICS GENERATED 
AXKXEKOEND) OF COMPLLATICN 4+¥44%* 


COMPILER STATISTICS: ELAPSED TIME 2043 SEC. 
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MTS FORTRAN H (OS REL 21.7) (C7204) 


COMPILER OPTIONS - NAME= MAIN, OPT=01,LINECNT=59, 
SOURCE, EBCDIC, NOLIST, DECK, NCLOAD, NOMAP 


G 
ISN 0002 SUBROUTINE INTRAN (E,N,Q) 
G INVERT AND TRANSPOSE MATRIX E, OR TURN E UP 
Cc SIDE DOWN 
ISN 0003 COMPLEX E(N,N),Q(N,N) 
ISN 0004 Hoe50eis=15N 
ISN 0005 DOESOR I= 
ISN 0006 50 O(J,d3) = E«N+1-0,33) J 
ISN 0007 RETURN 
ISN 0008 END 
*OPTIONS IN EFFECT* NAME= MAIN, OPT=01,LINECNT=59, 
*OPTICNS IN EEFECT* SOURCE, EBCDIC, NOLIST, DECK, NOLOAD, NOMAP 
* STATISTICS* SOURCE STATEMENTS = 7 ,PROGRAM SIZE = 562 


Pot o te Coree NO DIAGNOSTICS GENERATED 
FSS ee OF CCONPILAT LON 24s es 


COMPIUENSSTATISTUCSs ELAPSED Tims si? LOMSEC. 
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MTS FORTRAN H (OS REL 21.7) (CT204) 


COMPILER OPTIONS - NAME= MAIN, OPT=01,LINECNT=59, 
SOURCE, EBCDIC, NOLIST, DECK, NOLOAD, NOMAP 


ISN 0002 SUBROUTINE INVV(AA,N,AZNV) 
Cc A ROUTINE TO INVERT VANDERMONDE MATRIX 

ISN 0003 COMPLEX A,AINV,AA,W 

ISN 0004 COMPLEX CMPLX 

ISN 0005 DIMENSION AB(N,N) ,AINV(N,N) ,A(20,40),ID (20) 

ISN 0006 NN=N+1 

ISN 0007 N2=24N 

ISN 0008 DOe100e1=1,-N 

ISN 0009 ID (1) =I 

ISN 0010 iG {MOO cHEST 4 

ISN 0011 100) AM aye (2,5) 

ISN 0012 DO 200 I=1,N 

ISN 0013 DO 200 J=NN, N2 

ISN 0014 MO) Gt aye 

ISN 0015 DO 300 I=1,N 

ISN 0016 SO0mA(r Nt) =1. 

ISN 0017 Kaa 

ISN 0018 1 CONTINUE 

ISN 0019 Dar (CABS (L (hy K))) 3, 99978 

ISN 0020 CuK he 

ISN 0021 DO 4 J=KK,N2 

ISN 0022 A(K, J) =A (Kd) 7a (KK) 

ISN 0023 Donde 16 N 

ISN 0024 IF (K-12) 41,4, 41 

ISN 0025 44 W=A (I,K) *A(K,J) 

ISN 0026 A (I,J) =A (I,d)-® 

ISN 0027 IF (CABS(A(ZI,dJ))-.0001*CABS(W)) 42,4,4 

ISN 0028 OZ eeAa ee CME L (Oe Cl) 

ISN 0029 4 CONTINUE 

ISN 0030 K=KK 

ISN 0031 IF (K-N) 1,2,5 

ISN 0032 Sa) ek Al ait 

ISN 0033 DOM 17 

ISN 0034 IF (ID (J) -I)11,8,11 

TSN 0035 8 DO 10° K=1)N 

ISN 0036 AINV (I,K) =A(d,N+K) 

ISN 0037 10 CONTINUE 

ISN 0038 19 CONTINUE 

ISN 0039 12 CONTINUE 

ISN 0040 RETURN 

ISN 0041 999 PRINT 1000 

ISN 0042 41000 FORMAT(19H MATRIX IS SINGULAR) 

ISN 0043 RETURN 


ISN 0044 END 


= aye ae 





FOPEPEONS ENEEPEECT= NAME= MAIN, OPT=01,LINECNT=59, 
SOR DLONOm UNE Pr h GL SOURCE, EBCDIC, NOLIST, DECK, NOLOAD, NOMAP 
Socal top eco* SOURG lot fo CENTS = 43 ,PROGRAM SIZE = 


FOLATISTICS*=) NOP DESGNOSTICS GENERATED 


444k END OF COMPILATICN ¥%%%*% 


COMP TUPR > LA RS neo. ELAPSED TINE Fey CNW C7 
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Gao Nesnx«ecutvon Program 
This program reads in the data from file Data and 
solves the equations using the fast algorithm in (9.5.5). 


The whole program has only two control cards: 


SRUN COMPUT 7=DATA 


SENDFILE 


Output of the program is listed in next page. The 
result shows that the total execution time in solving a 
linear system of 64 complex equations (128 unknowns) has 
been reduced to 2.18 seconds. 

The maximum relative error occurs in the imaginary 


value of 62 where 


10000000" = 0. 9999946 


1.0000000 = 0.00054% 


Max Rel Error= 


Toa 
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¢ open deer dt basen Rea 2H ao daqtu0 


. paiclea at otis ealotboe* led eds See eee Jloeet 
7 7 
“Ad (eswoeket BET) “ses ae ee he 2o: modem naeakt 


72 :  labneeee Si, 3 43 Beoubes. need: 


wanes « ‘fa 
ft id Hr a2Yous ovens afeyz ptiine’ sks” oat ae 


: uf Di 4s << « Pee es a { va 
PEQUG Qo rt eee = San ans et 7 


ThemVarlUeSa © feet cr ct 


PE = Fy) The expected values should be either 241 


or O+ OF 


Fs (eel = 
Rat Ej 
E32) = 
f(a) = 
F( 5)= 
EGS i= 
1 DV 
F( 8) = 
Wet 3) = 
F(10) = 
F(11)= 
PGEl2— 
GES) = 
F (14) = 
F(15)= 
F (16) = 
EQ} 
F (18) = 


¥F (19) = 


di 


0.0000018 
0.0000026 
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APPENDIX IV 


DIGITIZATION OF A SCATTERING FORMULA 


The theoretical derivation in Chapter VII arrived 


at the formula (7.1.7), which is 


HO) = i. G(X"/X) £,(X) ul (x) av (A4.1) 


where Uo, (X") iseines optical) fieldsscattered by a weakly 
scattering object of microscopic size 
£(%) Het hemscatrLeningepolLentiale ot thesobyect 
uy (X) Tse aneincvacnt eeteld 


G(X0/xX weisean Green se sunctromechosen sto be 





~exp (ik, |X'-X|)/4m k'-X 
Integration over the volume V will first be per- 
formed inside each blocklet and then be summed up block- 
tetyby blockler: 


eae) = L . G(Xl 7x) fo (x) uo (x) dv . (Adm) 


Putting in the far-field approximation Givens nee eeles.) 7 
LiCl ilar ee LO) 


-exp (ik, |X']) 





en Aa) = } L. exp (-ik+X) u, (X) £ , (X) av. 
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(Ades) 


Let the vector X be resolved into two components 


Xe and sox. 
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bilochle frase mincdetine (ie2.1)) moxswillabera Vector 
Grommet hemcen ero met nea tne block ler pointing CO, any, 


HOLnEeXeWwi eh nethnemven blockler. 


sgt 6") 
eG Sa | exp [-ik- (X_+6X) ] 
ah Am | Xan v ?V ¥ 
Vv 
Ua (X+6X) £, (X+6X) dV is. (A4.5) 


The incident Field uy (X +X) Wa lie ves split up 
i cOmeWOM Da Tass 
Oe (Xe 0X) 
Onn: 


a (X_+0X) = Lo (X_) ea = aoe (x) Os (6X) a (A4.6) 


Physical meaning of (A4.6) can be seen ieee Garis balCces, 


uy (X) is a plane wave propagated along the direct onnoL 


Wms Gia — igh Teleweeh CEIEIS 
uy (%) = exp (ik) 2) ; (A4.7) 
U(X +0X) = exp [ik, (z2,+92) ] ; (A4.8) 
Toh TOs Bap )ALL (0.5) St ickgs nh foray) FEF ESPs < (Aa, 9) 


uy (X+5%) may be considered as separable into two waves: 


uy (X) acts like a carrier wave, which is actually the 


incident wave uy (X) measured at the center of the vth 


blocklet. wu, (5X) is a modulator wave, which looks after 
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the phase difference within the vth blocklet. It shows 
that in this case u, (6X) is independent of v. 

From (A4.5), and writing £ (x +6X) as £, (6X) since 
x, LSsaconstancerorethe integral over Net 


-exp (ik, |[X"]) 


BRE) ) exp (-ik-X_) u)(X,) | 


amex Ul Vv V 
Vv 


exp (-iK- 6X) uy (6X) £, (6X) dav. (A4.10) 


Define toy as a weighted statistical average of 


f, over the vth blocklet: 


Jv, Shen CO aces wee dv 


fe (X08 ec . (A4.11) 


| exp (~iK. 6X) up, (5X) adv 
V 
Vv 


In many cases, i (%) varies very slowly within the tiny 


space of a blocklet. We can then assume 


pooner) ee (A4,12) 


Then the left hand side of (A4.11) becomes £i(%)- The 
numerator on the right hand side of (A4.11) is the 
integral ing (A410) eLetwusecal thes denominater oF 
CASI CCK) 


a(K) = | exp (~iK-6X) u, (SX) ay. (A4.13) 
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When the incident wave Uo is a plane wave along 


the Z-direction, 
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o(K) AV, where $(K) is called a form factor. 
(A4.14) 


This shows that $9(K) is independent of v. It is the 
same for all blocklets. (A4.14) also shows that when 
Ax, Ay, AZ approach zero, i1.¢€. the blocklets become 


very small, $(K) approaches 1, a(K) approaches AV. 9(K) 
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Since ¢(K) and AV are both independent of v, they can 


be taken out of the summation symbol. 
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Let 
eee Xie | 
F(K) = ——=°-—____——-__ ,_ K = k, Kee (AAT) 
exp (ik, |X" |) $(K) AV eh 
Beh = ee) OC Bees (A4.18) 
Finally, (A4.16) becomes 
ip () =) ERS) ase) (A4.19) 
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APPENDIX Y 
PROPAGATION OF LIGHT IN A QUASI-HOMOGENEOUS *.DIUM 
Letetnesbenavioureon Jiqht in a quasi-homogeneous 


medium be described by the Maxwell equations (A5.1) and 


Material equations (45.2) sands (A5.3) 


3 Sate. 
curl H - = ae D = Er Jor 
Tenet a 
CuGnLr ES a ae Be=-100, (Asis) 


jo) eS as (A5.2) 


Bia lier. (A5ne) 


For a non-conductive and non-magnetic medium, the elec- 
tric charge density op and the electric current density 
jmare bochezero, and the permeabilaty wis equal eto 


Ue Vee nce 


curl H - = < De=2 Or, (A5. 4a) 
Lec _ 

Cup lees, G At H = O ’ (A5.4b) 

dive tee) =AOV: (A5.4c) 


Applying the Operator CULL, CO mbOthe Sides sor 


(A5.4b), 
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eibteih velehell Aeo ss 4 Abhel Wl Sa ser - (A5.5) 


PU Cl igen hoe aa meteor (A> ..0):, 
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Subs eileutange (socal nco (A>. 0); 
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ybbelk Tekbeell fo) Ge =; meee 


Fe Ome (A5.7) 
Cio 


Remember the mathematical relations between the 
vector operators Curl, div, and grad) (see, for cxample, 
Pipe: Mathematics for Engineers and Physicists. Chap.15, 


D594 =McGraweHEl wo 58)): 





CligmeCcuml mie mg racdecd i Vaiis— V-E (A5.8) 
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V-E - grad div E + —> B= 07. (A510) 
Ca Neo 
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Equation (A5.11) describes, in exact mathematical 
terms, the optical wave E propagated in a non-homogeneous, 


non-conductive, non-magnetic medium. Since the material 
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is non-homogeneous, ¢€« should be a function of (x,y,z). 
For some physical problems, ¢« varies very slowly in 
comparison with the variation of an optical field E 
in Space. in this case, the last term in (A5.11) is 
negligibly small and Mave ben drouped a ON Ustary eras 
SiMnpuiErCcet on me letmuceestimate thes ectual Valuessof 
the three terms in (A5.11). As a very rough Spey 
we shall assume that E is a plane wave. 

Section 7.la gives E = Refu exp(-iwt)]. In the 


case of a plane wave propagated in z-direction, 





E = Re[A exp (-iwt ~ik 2) ] ‘ ; (AS sZ) 
Then 
2 9°E o 2 
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meterse-, Orshasea Similar Magnitude as the first 
erin. 
Regarding the third term grad (= 6, Pepe vel Fs) ial 


(A5Se lee em cevarres slowly Unespace, we Can expect 


14) sve eregeinae => ie Vee 
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APPENDIX VI 


SOLUTION OF A NON-HOMOGENEOUS HELMHOLTZ PARTIAL 


Sc 


Ors 


DIFFERENTIAL EQUATION 


Wemare goingeco solve, the equation an. (/.1.6) for 


The equation takes the form: 


V-u + ae =) £5 (AG ay) 
sc Ons e SO 

on Re oe ete (A6.2) 
sc OSE S oO 


Since the scattered field is due to multiple point 


sources of density fouoe contained within the volume V, 


we shall choose the point-source Green's function 





expo key) Xe xa) 
Ce ee (A6. 3) 
Po 


This Green's function has two properties: 


il 


hey 


Differentiating (A6.3) gives 


dG 


i IS ey Te eee ae PS (A6.4) 
a|x'- x| 


THLSeGLeenusmLtunceLonsalmsuLes 


V-C 4+ ko G = 6(x!- x) eNe ae 
qe = MES jay Ke G. (A6 .6) 


205 










york ‘axetaiovate 


7 


: 8s 
> feel. a (jewpa <a Ww Is@ ag sk oe al on } 
incest ano ‘bed =a noasnaye sate 





a? 2 * ad Wr «ens Bint 3 ia 
et SL Soe 398 wd 20ak8 


Ci] LEB OD", « re tolemeb ie icone 
tl Hse separa hee, Seis papeds: mero = : a Bak: 


tap Be ; 4 UW Gaéen not tsaut eases: 
oar se Pe LORS RRR TD 





We shall make use of the second form of Green's 
theorem, which states 
| (Gy-u -u y°G) av = (GVYu__- u__VyG)-dA (A6.7) 
Sc sc = sc eye! : 2 
V 
Them eetenandestde Of (AG), On ssubStilution 
Wi CMe Oreo amano A6.0)), O1lvVes 
2 
| [-Gk~u sriGhe Abl Se Abi 
Ogsc sO 
V 


2 { 
eye toe us 6(% Oehavi 


# ac fujlav- uj) . (AG .8) 


The rignt hand side of —(AG.7), taking the gradients 


Tne herOUEeCct Lone Olms(X. —X), -DeComes 


€ et ee es ee (A6.9) 
ae) sc 9(X'-X) Ane 


Substitution of (A6.4) turns (A6.9) to 


egies ‘vee “se 1 laa (AGEEO} 
7 (KPAX) We Paice \ : ° 


We have the radiation condition (see, for example, 


BuLckovelJ6S mp swoorand p.ol7)} 


du 
pei eparesep eee ete a Oya GL) |X'=X|>0. (A6.11) 
Also, 
Ysc 
rad —-oO, Se hk x er co ames (A6.12) 


Theretore (AGvL0), or the right hand side ols (AG.7/)), 


disappears under far-field condition (Gia lena yes 


206 













— te 


‘aes 40 ia Sache bd ae. 


aA} . ADC oa “ewe 
noite tvecdie ao 7 + ae is ebfte Soa ved ey 7 
ugh, (azaA) Bh het Pie 


. “| z _ : : 
vil (ks 200. a ed ee a 28° oar ood | "ge 
: 2 -o De ? : 


- 





4 hee a 
mek av a a 


‘ 
<a 


ld perset - [t ao ohite brine S.0gks fT! PS ee 


sqmaod” , (eet, “2g re snd nk i. 


as 
ip 


4m} 
a 
(a, spa) )AGky aes a ae — 
\\ fae pitr niteg peer rotgpe = (CS 


Soon TE uted oe | oe aul ree 


' ar ee 7 


267 


Equating the left hand Pen eee) LOEZEXLO;, 


we have 
u = | te aw Gohve CA GeeloD 
sc Se) 


Witchmeicmi hems requined solution ( /l./) (oO equation 
(FJeleo)y. This solution expresses the scattered field 
ue in terms of the source distribution FY, for any 


pointe Xie insidesor outside the scatterer, in particular 


Couecine the tar-nield) region. 
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